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NUMERICAL ANALYSIS OF A ROBUST FREE ENERGY 
DIMINISHING FINITE VOLUME SCHEME FOR PARABOLIC 
EQUATIONS WITH GRADIENT STRUCTURE 

CLEMENT CANGES AND CINDY GUICHARD 


Abstract. We present a numerical method for approximating the solutions 
of degenerate parabolic equations with a formal gradient flow structure. The 
numerical method we propose preserves at the discrete level the formal gradient 
flow structure, allowing the use of some nonlinear test functions in the analysis. 
The existence of a solution to and the convergence of the scheme are proved 
under very general assumptions on the continuous problem (nonlinearities, 
anisotropy, heterogeneity) and on the mesh. Moreover, we provide numerical 
evidences of the efficiency and of the robustness of our approach. 


1. Introduction 

Many problems coming from physics (like e.g. porous media flows modeling [11, 
10, 35]) or biology (like e.g. chemotaxis modeling [67]) lead to degenerate parabolic 
equations or systems. Many of these models can be interpreted as gradient flows in 
appropriate geometries. For instance, such variational structures were depicted for 
porous media flows [86, 70, 32], chemotaxis processes in biology [19], superconduc¬ 
tivity [5, 4], or semiconductor devices modeling [84, 68] (this list is far from being 
complete). 

Designing accurate numerical schemes for approximating their solutions is there¬ 
fore a major issue. In the case of porous media flow models — used e.g. in oil¬ 
engineering, water resources management or nuclear waste repository management 
— the problems may moreover be highly anisotropic and heterogeneous. As an 
additional difficulty, the meshes are often prescribed by geological data, yielding 
non-conformal grids made of elements of various shapes. This situation can also be 
encountered in mesh adaptation procedures. Hence, the robustness of the method 
w.r.t. anisotropy and to the grid is an important quality criterion for a numerical 
method in view of practical applications. 

In this contribution, we focus on the numerical approximation of a single non¬ 
linear Fokker-Planck equation. Since it contains crucial difficulties arising in the 
applications, namely degeneracy and possibly strong anisotropy, the discretization 
of this nonlinear Fokker-Planck equation appears to be a keystone for the approxi¬ 
mation of more complex problems. 

1.1. Presentation of the continuous problem. Let H be a polyhedral con¬ 
nected open bounded subset of (d = 2 or 3), and let tf > 0 be a finite time 
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horizon. In this contribution, we focus on the discretization of the modei probiem 


( dtu - V- (? 7 (u)AV(p(m) + I/)) = 0 in Qt,. := O x 
(1) < {ri{u)A'V{p{u)+ V)) ■ n = 0 on9nx(0,tf), 

Uu=o=“o infl, 


which appears to be a keystone before discretizing more compiex probiems. We do 
the foiiowing assumptions on the data of the continuous probiem ( 1 ). 

(Al) The function rj : IR+ —)■ IR_|_ is a continuous function such that 77 ( 0 ) = 0, 
r]{u) > 0 if u 7 ^ 0 and 77 is non-decreasing on K_|_. The function 77 is 
continuousiy extended on the whoie K into an even function. It is called 
the mobility function in reference to the porous media flow context. 

(A2) The so-called (entropy) pressure function p G is absolutely contin¬ 

uous and increasing on ( 0 , -boo) (i.e., 0 < p' £ iioc(( 0 i +oo)))i satisfies 
lim„_j.+oo p(w) = -1-00. In the case where p(0) = lim„^oP(w) is finite, the 
function p is extended into an increasing absolutely continuous function 
p : K — 7 ^ M defined by 

(2) p{u) = 2p(0) — p{—u), Vm < 0. 

We denote by 

if p( 0 ) = —00, 
if p( 0 ) > —00. 

and by Ip its closure in K. We additionally require that the function 
yljp' belongs to (and is in particular integrable near 0) and that 

lim„^o Vv{u)p{u) = 0. 

(A3) The tensor field A : 11 —> (L°°(K))‘^^‘^ is such that A(a;) is symmetric for 
almost every a; G 12. Moreover, we assume that there exist A* > 0 and 
A* G [A*, -boo) such that 

(3) ^*|vp < A(a:)v • v < A*|vp, Vv G for a.e. a: G 12. 



A is called the intrinsic permeability tensor field still in reference to the 
porous media flow context. 


(A4) The initial data ug is assumed to belong to L^(fl). Moreover, defining the 
convex function T : Ip —)■ IR+ (called entropy function in the following) by 

(4) T{u) = (j){a)-p{l))da, Vtt G Ip, 


we assume that the following positivity and finite entropy conditions are 
fulfilled: 


( 5 ) 


uq > 0 a.e. in 12 , 


uoda; > 0 , 


r(uo)da; < -boo. 


(A5) The exterior potential V : n —>■ K is Lipschitz continuous. 
Throughout this paper, we adopt the convention 
( 6 ) r(M) = -boo if p( 0 ) = —00 and u < 0 . 
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In order to give a proper mathematical sense to the solution of (1), we need to 
introduce the function ^ : IR_|_ —>■ IR_|_ defined by 

l>u _ 

(7) ^{u) = / A/? 7 (a)p'(a)da, Vu > 0. 

Jo 

Note that ^ is well defined since we assumed that belongs to LlM- More¬ 
over, in the case where p(0) is finite, then the formula (7) can be extended to the 
whole M, leading to an odd function. We additionally assume that the following 
relations between 77 and F hold: 

(A 6 ) There exists C > 0 such that 

( 8 ) 0 < ^(u) < C (1-I-r(u)), Vu G [0,-boo). 

Moreover, we assume that 

, . r(u) 

(9) —^ ^ -boo as u —^ -boo, 

r]{u) 

and that the function 

( 10 ) ^/^ ^ ^ is uniformly continuous on the range of 

Definition 1 (weak solution). A measurable function u is said to be a weak solution 
to problem (1) if 

i. the functions u and rj^u) belong to L°°(( 0 , tf); L^(f 2 )); 

ii. the function f{u) belongs to L‘^{{0,tf)',H^{^)); 

iii. for all function tp G x [0,tf);K), one has 



udttjj dxdt -b 


Uo’ipi', 0 )da; 



?7(u)AVy • ^tpdxdt 


AV^(u) • \/r]{u)'V'tjjdxdt = 0. 


Following the seminal work of [2], there exists at least one weak solution u to 
the problem (1). Denoting by 


PU 

(j){u) = / ri{a)p'{a)da, Vu G Ip, 

Jo 

the uniqueness of the solution (and even a L^-contraction principle) is ensured as 
soon as r] o (j)~^ £ ( 70 . 1/2 g^g g^jg^ j-yj ^ giightly weaker condition in the 

case of a smooth domain D). Moreover, u belongs to C([0, if]; I^(f^)) (cf. [31]) and 
u{-,t) > 0 for all t G [0,tf] thanks to classical monotonicity arguments. 


Remark 1 . 1 . Assumptions (K\)-(Ado) formulated above deserve some comments. 

• First of all, let us stress that Assumptions (Al)- (A2) and (A6) are sat¬ 
isfied if r]{u) = u and p{u) = log(u) as in the seminal paper of Jordan, 
Kinderlehrer, and Otto [65]. One can also deal with power like pressure 
functions p(u) = , but only for m > 1. Our study does not cover the 

case of the fast-diffusion equation m < 1 with linear mobility function (see 
e.g. [ 86 ]j because of the technical assumption (A6). 
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• The most classical choice for the mobility function 77 : IR_|_ —>■ K+ is r](u) = 
u. In this case, the convection is linear. In this situation, the formal 
gradient flow structure highlighted in %1.2 can be made rigorous following 
the program proposed in [65, 86 , 1, 3, 77] and many others. The gradient 
flow structure can also be made rigorous in the case where t] is concave 
(cf. [43] j and in the non-degenerate case r]{u) > a > 0. 

• One assumes in (Al) that rj is nondecreasing on K+. This assumption is 
natural in all the applications we have in mind. However, it is not manda¬ 
tory in the proof and can be easily relaxed: it would have been sufficient to 
assume that there exists 7 > 0 such that 

—^ vis), V[a,5] C K. 

This relation is clearly satisfied with 7 = 1/2 when rj is nondecreasing. 

• In Assumption (A2), the condition liin„_>oo p(m) = +00 ensures that 

r(M) 

lim - = +c», 

u—>-oo 'll 

whereT was defined in (AA). Given a sequence (««)„ C L^{ri) with bounded 
entropy, i.e., such that is bounded, then (Mn)„ is uniformly 

equi-integrable thanks to the de La Vallee Poussin’s theorem [42]. There¬ 
fore, a sequence (ura)„ with bounded entropy relatively compact for the weak 
topology of {LI). 

• Since the unique weak solution to the problem remains non-negative, the 
extension of rj and p on the whole K could seem to be useless. However, 
in the case where J3(0) is finite, the non-negativity of the solution may not 
be preserved by the numerical method we propose. The extension of the 
functions rj and p on K_ is then necessary. 

• Only the regularity of the potential V is prescribed by (Ab). Confining 

potentials like e.g. V{x) = for some a;* € fl, or gravitational 

potential V{x) = —g ■ x, where g is the (downward) gravity vector can be 
considered. 


1.2. Formal gradient flow structure of the continuous problem. Let us 

highlight the (formal) gradient flow structure of the system (1). Following the path 
of [86, §1.3] (see also [84, 87]), the calculations carried out in this section are formal. 
They can be made rigorous under the non-degeneracy assumption 77(11) > a > 0 for 
all It > 0. 

Define the affine space 


m = 


|it : D —K 


u{x)d.x 


/ uo{x)dx 

Jn 


of the admissible states, called state space. 

In order to define a Riemannian geometry on 911, we need to introduce the 
tangent space TfUfl, given by 


T„91l =\w:VL 


/ w{x)dx = 0 

Jn 
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We also need to define the metric tensor x —) 

a scalar product on T^Ai (depending on the state u) 

(12) 9u{wi,W2) = / cj)iW2dx= / wi^2da; = / /^(m) V(/)i • AV(()2 dx, 

Jo, Jo Jo 

for all wi,W 2 € 'Til®!, where (j)i are defined via the elliptic problem 

—'V-(ji{u)A'V(l)i) = Wi in n, 
r]{u)A'V(j)i • n = 0 on dfl, 


which consists in 


(13) 


I dx = 0. 


Note that T„S[Jl does not depend on u (at least in the non-degenerate case), but the 
metric tensor 0«(', •) does. So we are not in a Hilbertian framework. 

Define the free energy functional (cf. [64]) 


(14) (£ : 

and the hydrostatic pressure function 

f): 


—?► M U {+oo} 

u !->■ £(u) = / (r(u(x))-I-u(x)y(x)) dx, 

Jfi 


Ip X n 


(u, x) i-> \)(u,x) = p(u) -I- V[x) = Du<B{u). 

Then given w € TyfXft, one has 

(15) Du€{u) ■ w = / f)(u(x), x)?ii(x) dx = / ?7(u)V[)(u, •) • AV0 dx, 

Jq Jq 

where (j) is deduced from w using the elliptic problem (13). Moreover, thanks to (12), 
one has 


(16) 


Qu{.dtu,w) = / dtucfdx, \/w e 
JQ. 


m. 


In view of (15)-(16), the problem (1) is equivalent to 

(17) Qu{dtu, w) = -Du€{u) ■ w = -0„(Vui£(u), w), Vw G I„®1, 

where the cotangent vector Du^{u) G (TuSJt)* has been identified to the tangent 
vector Vii£(m) G TuWl thanks to Riesz theorem applied on T^Wl with the scalar 
product 0u. This relation can be rewritten as 

(18) dtu =-Vu£(u) = V-(p(u)AVl)(u, ■)) in TuM, 

justifying the gradient flow denomination. 

Choosing w = dtu in (17) and using (18), we get that 

^€(u) = Du€{u) ■ dtU = -Dy,e{u) ■ V„€(m) = [ r]{u)Vi){u, •) • AV[)(m, ^dx. 
ac Jq 

An integration w.r.t. time yields the classical energy/dissipation relation: Vt G 

[0,tf], 

(19) £(«(•, t)) - (£(uo) 

+ / / i?('u(x, t))A(x)V()(u(x, r), x) • V(](u(x, r), x)dxdT = 0. 

Jo Jq 
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The fact that a physical problem has a gradient flow structure provides some 
informations concerning its evolution. The physical system aims at decreasing its 
free energy as fast a possible. As highlighted by (19), the whole energy decay 
corresponds to the dissipation. As a byproduct, the free energy is a Liapunov 
functional and the total dissipation (integrated w.r.t. time) is bounded by the free 
energy associated to the initial data. The variational structure was exploited for 
instance in [86, 21, 22, 93] to study the long-time asymptotic of the system. 

1.3. Goal and positioning of the paper. The goal of this paper is to propose 
and analyse a numerical scheme that mimics at the discrete level the gradient flow 
structure highlighted in §1.2. Since the point of view adopted in our presentation 
concerning the gradient flow structure is formal, the rigorous numerical analysis of 
the scheme will rather rely on the well established theory of weak solutions in the 
sense of Definition 1. But as a byproduct of the formal gradient flow structure, the 
discrete free energy will decrease along time, yielding the non-linear stability of the 
scheme. 

There are some existing numerical methods based on Eulerian coordinates (as 
the one proposed in this paper). This is for instance the case of monotone dis¬ 
cretizations, that can be reinterpreted as Markov chains [81, 45], for which one 
can even prove a rigorous gradient flow structure. Classical ways to construct 
monotone discretizations in the isotropic setting A(a:) = X{x)ld are to use finite 
volumes schemes with two-points flux approximation (TPFA, see e.g. [51, 50]) or 
finite elements on Delaunay’s meshes (see e.g. [37]). An advanced second order 
in space finite volume method was proposed in [18] and discontinuous Galerkin 
schemes in [80, 79]. However, these approaches — as well as the finite difference 
scheme proposed in [78] — require strong assumptions on the mesh. Moreover, 
the extension to the anisotropic framework of TPFA finite volume scheme fails for 
consistency reasons (cf. [50]), while finite elements are no longer monotone on a 
prescribed mesh for general anisotropy tensors A. In [61, 56], it appears that all 
the linear finite volume schemes (i.e., schemes leading to linear systems when linear 
equations are approximated) able to handle general grids and anisotropic tensors 
may loose at the discrete level the monotonicity of the continuous problem. 

The monotonicity at the discrete level can be restored thanks to nonlinear cor¬ 
rections [27, 71, 30, 72]. Another approach consists in designing directly monotone 
nonlinear schemes (see, e.g., [66, 92, 75, 76, 44, 90]). However, the monotonicity of 
the method is not sufficient to ensure the decay of non-quadratic energies. More¬ 
over, the available convergence proofs [44, 30] require some numerical assumptions 
involving the numerical solution itself. Finally, let us mention here the recent con¬ 
tribution [59] where linear monotone schemes are constructed on cartesian grids for 
possibly anisotropic tensors A. 

In the case where A = and r]{u) = u, the formal gradient flow structure can 
be made rigorous in the metric space 



endowed with the Wasserstein metric with quadratic cost function. Several ap¬ 
proach were proposed in the last years for solving the JKO minimization scheme 
(cf. [65, 3]). This requires the computation of Wasserstein distances. If d = 1, 
switching to Lagrangian coordinates is a natural choice that has been exploited for 
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example in [69, 20, 82, 83]. The case d > 2 is more intricate. Methods based on a so- 
called entropic regularization [13, 88] of the transport plan appear to be costly, but 
very tractable. Another approach consists in solving the so-called Monge-Ampere 
equation in order to compute the optimal transport plan [15]. Let us finally men¬ 
tion the application [14] to the Wasserstein gradient flows of the CFD relaxation 
approach of Benamou and Brenier [12] to solve the Monge-Kantorovich problem. 

Motivated by applications in the context of complex porous media flows where 
irregular grids are often prescribed by the geology, the scheme we propose was 
designed to be able to handle highly anisotropic and heterogeneous diffusion tensor 
A and very general grids (non-conformal grids, cells of various shapes). It relies 
on the recently developed Vertex Approximate Gradient (VAG) method [52, 54, 53, 
26], but alternative versions can be inspired from most of the existing symmetric 
coercive methods for approximating the solutions of linear elliptic equations [61, 56]. 
Moreover, we want our scheme to mimic at the discrete level the gradient flow 
structure highlighted in §1.2. This ensures in particular the decay of the discrete 
counterpart of the free energy, and thus the nonlinear stability of the scheme. 

A nonlinearly stable control volume finite elements (CVFE) scheme was proposed 
in our previous contribution [34] (see also [33]). The nonlinear CVFE scheme [34] is 
based on a suitable upwinding of the mobility. It is only first order accurate in space 
while linear CVFE schemes are second order accurate in space. Moreover, it appears 
in the numerical simulations presented in [34] that this nonlinear CVFE scheme 
lacks robustness with respect to anisotropy: its convergence is slow in particularly 
unfavorable situations. 

The main goal of this paper is to propose a scheme that preserves some impor¬ 
tant features of the one introduced in [34] (possible use of some prescribed nonlin¬ 
ear test function, decay of the physically motivated energy, convergence proof for 
discretization parameters tending to 0), without jeopardizing the accuracy of the 
scheme compared to the more classical approach based on formulations with Kirch- 
hoff transforms (see for instance [57, 89, 58, 9]). Convincing numerical results are 
provided in §5 as an evidence of the efficiency of our approach. Two theorems are 
stated in §2.4 (and proved in §3 and §4) in order to ensure the following properties: 

(1) Theorem 2.3. At a fixed mesh, the scheme, that consists in a nonlinear 
system, admits (at least) one solution. This allows in particular to speak 
about the discrete solution provided by the scheme. Moreover, we take ad¬ 
vantage of the gradient structure of the scheme for deriving some nonlinear 
stability estimates. 

(2) Theorem 2.4. Letting the discretization parameters tend to 0 (while con¬ 
trolling some regularity factors related to the discretization), the discrete 
solution converges in some appropriate sense towards the unique weak so¬ 
lution to the problem (1) in the sense of Definition 1. 

Remark 1.2. We only consider potential convection in the paper. A more general 
convection with speed v can be split into a potential part and a divergence free part: 

V = —W -I- V X w. 

We suggest for instance to use classical (entropy stable) finite volume schemes (see 
for instance [74]^ for the divergence free part combined with our method for the 
potential part. 
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2. Definition of the scheme and main results 

As already mentioned, the scheme we propose is based on the so-called VAG 
scheme [52], In §2.1, we state our assumptions on the spatial mesh and the time 
discretization of (0, U). Then in §2.2, we define the nonlinear scheme we will study 
in this paper. The gradient flow structure of the discretized problem is highlighted 
in §2.3, where a variational interpretation is given to the scheme. Finally, in §2.4, 
we state the existence of discrete solutions to the scheme and their convergence 
towards the unique weak solution as the discretization parameters tend to 0. 

2.1. Discretization of Qt^ and discrete functional spaces. 

2.1.1. Spatial discretization and discrete reconstruction operators. Following [53, 
26], we consider generalized polyhedral discretizations of 11. Let M be the set of 
the cells, that are disjoint polyhedral open subsets of 11 such that UkgA 4 ^ 

Each cell k G At is assumed to be star-shaped with respect to its so-called center, 
denoted by a:^. We denote by J- the set of the faces of the mesh, which are not 
assumed to be planar if d = 3 (whence the term “generalized polyhedral”). We 
denote by V the set of the vertices of the mesh. We denote hy Xg G id the location 
of the vertex s G V. The sets Vk, and denote respectively the vertices and 
faces of a cell k, and the vertices of a face a. For any face a one has Vo- C Vk- 

Let Afs denote the set of the cells sharing the vertex s. The set of edges of the 
mesh (defined only if d = 3) is denoted by S and 8a denotes the set of edges of the 
face a G J-, while denotes the set of the edges of the cell k. The set Ve denotes 
the pair of vertices at the extremities of the edge e G 8. In the 3-dimensional case, 
it is assumed that for each face a G J-, there exists a so-called “center” of the face 
Xa such that 

(20) Xa — ^ ] /3o-,s^S 5 with ^ ] P( 7 ,s ~ 1; 

sGVct sGVct 

and Pa.s > 0 for all s G Vo-- The face cr is then assumed to match with the union of 
the triangles Ta^e defined as by the face center Xa and each of its edge e G 8a. A 
two-dimensional example of mesh M is drawn on figure 1. 

The previous discretization is denoted by V, and we dehne the discrete space 

Wry = {v = e . 

In the 3-dimensional case, we introduce for all cr G A the operator la : Wt> -G M 
defined by 

Itriv) = ^ l3a,sVs, Vt) G Wp, 

sGV,, 

yielding a second order interpolation at Xa thanks to the dehnition (20) of Xa. 

We introduce the simplicial submesh T (a two-dimensional illustration is pro¬ 
vided on Figure 2) defined by 

• T = {Tt^,a, K G Al, cr G Ak} in the two-dimensional case, where T^^a denotes 
the triangle whose vertices are a;„ and Xs for s G Vo-; 

• T = G M.,a G As, e G 8a} in the three-dimensional case, where 

Tfi.a.e denotes the tetrahedron whose vertices are Xa and Xg for s G Vg. 
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Figure 1 . The primal mesh A4 can be made of cells with various 
and general shapes. Degrees of freedom are located at the so-called 
cell center (blue dots) and at the so-called vertices (red dots). 



Figure 2. The simplicial submesh T is derived from the primal 
mesh A4 by decomposing the primal cells k G M into triangles if 
d = 2 or tetrahedra if d = 3. This construction is possible since k 
was assumed to be star-shaped with respect to a ball centered in 
Xk- 


We define the regularity Oj- of the simplicial mesh T by 

(21) d7- = max—, 

tgT pt 

where and px respectively denote the diameter of T and the insphere diameter 
of T. We denote by 

(22) hr = niaxhT 

tgT 

the maximum diameter of the simplicial mesh. We also define the quantities and 
£s quantifying the number of vertices of the cell n and the number of neighboring 
cells for the vertex s respectively: 

(23) 


4 = #V«, 4 = Vk gM,VsG V. 
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This allows to introduce the quantity 

(24) £x) = max ( max ) , 

\/3GA4UV ^ kGAI / 

controlling the regularity of the general discretization V of 12. 

Denoting by Hj- C the usual Pi-finite element space on the simplicial 

mesh T, we define the reconstruction operator ttj- : W-v —t Hj- by setting, for all 
V e Wxi and all (s, k, tr) G V x A4 x J", 

(25) Trrv{xs) = Vs, tttv{xk) = and 7Trv{x^) = Idv). 

This allows to define the operator Vr : W-d —t (T°°(12))‘^ by 

(26) Vrii = Vttt-i), Vu G Wt>- 

We aim now to reconstruct piecewise constant functions. To this end, we intro¬ 
duce a so-called mass lumping mesh T) depending on additional parameters that 
appear to play an important role in practical applications [54]. Let k G A4, then 
introduce some weights (oK,s)sgy such that 

(27) > 0, and ^ <1, Vs G V„. 

sGV,^ 

Denoting by meas(K) = J], da; the volume of k, then we define the quantities 
= Q!„,smeas(K), VkGA4,VsGVk, 

m«; = meas(K) - I]sgv„ Vk G M, 

so that one has 

TO/3 = meas(12). 

/3eMuv 

For all K G A4, we denote by uJk. and some disjointed open subsets of k, such 
that 

U I IJ 

\sGV„ 

and such that 

meas(a;K) = to^ and meas(a;K_s) = 'mK,s, Vk G A4, Vs G Vk- 

Note that such a decomposition always exists thanks to (27)-(28). Then we denote 
by 

~ € V. 

kGMs 

The mass-lumping mesh T) is made of the cells and ojs for k G A4 and s G V. An 
illustration is presented in Figure 3. In [54], a study focused on the repartition of 
the porous volume between nodes and centers is proposed, in the case of a coupled 
problem (transport of a species in the flow of a fluid in a porous medium). The 
influence of this repartition mainly concerns the transport part, and not the fluid 
flow. In the framework of the present paper, we numerically observe that this 
repartition has not a strong influence in the cases studied. However, since the 
function rj can vanish, the limit choice to„ = 0 (resp. rrig = 0) can lead to singular 
cases, and therefore must be prevented. 
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Figure 3. The mass-lumping mesh V contains one cell and Ws 
by degree of fredom. The mass repartition depends on the factors 
Ok.s introduced at (27). 


In what follows, we denote by 


(29) 


Cv — 




/Sgaiuv TT-j-e^dx 


where e^,/3sAIUVis the unique element of W-d such that 


-Kr^pix-y) = SJ, G Mu V, 

and SJ is the kronecker symbol. 

We can now dehne the piecewise constant reconstruction operators nj) : Wt> 
L°° n BV(fl) and : Wx, ^ n BV(n) by 


(30) 

TTv{v)ix)= vMu^ix)+Y'^^ 

loTg (^) 1 

\/x G n, Vv G Wx>, 


k£A4 sGV 



and 




(31) 

tym{v){x) = Y 

Vx G 

Wv G W-D- 


iteM 


Let / : K —)■ K be a possibly nonlinear function, then denote by 

f{v) = «s)«;gA 1 .sgv e ■ 

Notice that in general, 

^rifiv)) -f- /(7rr(t’)) and Vr(/('u)) + V/(7rr(ii)), 

but that 

(32) 7r-D(/(v)) =/(7rx)(v)) and 7rAi(/(ii)) =/(7r^(f)), Vn G Wp. 

2.1.2. Time-and-space discretizations and discrete functions. Let TV > 1 and let 
0 = to < ^1 < • • • < tA-i < tA = tf be some subdivision of [0,tf]. We denote by 
At„ =tn — tn-i for all n G {1,... TV}, by At = (Ati,..., At^)^ G , and by 

(33) At = max At„. 

^ ' l<n<A 

The time and space discrete space is then defined by 

= {v = «,<),eA4.sGV,l<n<A ^ | . 
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For V e Wx),At and n € {1,..., N}, we denote by 

We deduce from the space reconstruction operators Tr-p, ttm and ttj- some time 
and space reconstructions operators TTp At, TT^ ^t, TTp Ai : Wp.At map¬ 

ping the elements of Wp^At into constant w.r.t. time functions defined by 

7rx).At'w(-,f) = 7rp(v”), and nT,A.tv{-,t) = irriv"^) 

iff S The gradient reconstruction operator Vp, At : Wp^At —>■ {L°°{Qtf))‘^ 

is then defined by 

^r,^tv = V'TTp.AtL’, Vp e Wp,Af- 


2.2. The nonlinear scheme for degenerate parabolic equations. For tt G At, 

we denote by A^, = (<is,s')s ^'ev ^ the symmetric positive definite matrix 


whose coefficients are defined by 

(34) <s'= / A(x)Vres(x) ■ Vres'(x)dx = 

J K 

It results from the relation 

■n']-e.^{x) + ^ ■n'fGsix) = 1, Va; G tt, Vk G At, 

that, for all u^v G Wp and all tc G At, one has 

(35) A{x)Vru{x) ■Vrv{x)dx = ^ - Us){vk - Vs')- 

sGV,,s'GV,, 

For K G M, we denote by : Wp —)■ the linear operator defined by 

(<5«;«)s =Vk- Vs, Vs G V«, Vp G Wp. 

With this notation, we obtain that (35) rewrites as 

/ A{x)V'ru{x) ■ V'j-v{x)dx = ■ A^S^u, Wu,v G Wp, Vk G At. 

J K 


In order to deal with the nonlinearities of the problem, we introduce the sets 
C Wp and CL Wt),A/ of the cLdmissiblG stcitcs defined by 

u G Wp‘^ iff Vi^Glp, Vi/GAtUV, 

and 

vew^% iff i;”GWg'^, Vn G {!,..., iV}, 
while we denote by W™ C Wp the set of finite entropy vectors: 

(36) u G IWp" iff £p(u) := / (r(7rpi;) + 7rp(i;)7rp(V)) da; < oo, 

an 

where V = (V^, t4)K,s G Wp is defined by 

(37) W« = W(x«,) Vs = V{xs) VK:GAt,VsGV. 

It is easy to check that, thanks to assumptions (A2) and to the definition (4) of 
the convex function F, one has Wp'^ C W™. 
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Given u G and V G Wxi, we define the discrete hydrostatic pressure h{u) = 
((!«;(««), G Wv by 

(38) f)K(u„) = p(u„) + 14, [}s(us) = p(ms) + V;, VKG7W,VseV. 


The initial data Uq is discretized into an element G by 

(39) itg = — [ uo{x)dx, \/f3 G MUV, 
so that, thanks to (5), 

(40) / TTx>{u°)dx = / uoda; > 0. 

Jq Jq 

Let us state a first lemma that ensure that the discretized initial data has a finite 
discrete entropy. 

Lemma 2.1. Let uq G L^(fl) be such that (A4) holds, V be such that (A5) holds, 
vP be defined by (39), and V defined by (37). Then there exists C depending only 
on ||uo||Li(n) and |!VF||ioo(n)d such that 

(41) 2:x>(m°) < 2:(uo) + Chj- < €{uq) + C'diam(n), 

where the entropy functional 6: is defined by (14) and its discrete counterpart €x) 
is defined by (36). In particular, vP belongs to Wffi. 

Proof. We deduce from Jensen inequality that 

r(u^) < — [ r(uo)da;, 

whence, thanks to (A4) and to the definition (36) of the discrete entropy func¬ 
tional C-pi one has 

£x)(m°) < £(uo) + / ue){'KT>'V-V)dxdx. 

Jo, Jo, 

The last term in the above inequality is equal to zero thanks to (40) . The Lipschitz 
regularity of V yields 

WmyV-Vh^^n) < ||Vy|U=o(a).hr, 

so that 

^v{u°) < C:(uo) -I- ||wo||Li(n)ll 

□ 


With all this setting, we can present the scheme we will analyze in this contri¬ 
bution. For u G P^S'^A.ty '"^0 introduce the notation 


(42) 




v{0+v{Us) 

2 


Vn G M, Vs G VK,Vn G {1,... ,N}. 


Given ii” ^ G IF™, the vector m” G is obtained by solving the following 

nonlinear system: 


u" — u 


n—1 




AG 


+ F,,s(w") = 0, 

sGV,. 


(43a) 


VkGM, 
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(43b) 


Atr, 


Y, = vsev, 


K.&Mb 


(43c) F«,s(m”) = Y “ f)s'«')), Vk € M, Vs € V«, 


s'GV» 


(43d) 


i^..s(«”) + i^s,.(«”) = 0, VAt e M, Vs G V«. 


The scheme (43) can be interpreted as a finite volume scheme, the conservation 
being established on the cells and ujg for k G A4 and s S V. As a direct 
consequence of the conservativity of the scheme, one has 


(44) 


f TTxi{u^)dx = j tt'd{u'^ ^)da;. 

Jfi Jq 


However, contrarily to usual finite volume schemes, the fluxes are not issued 
from the computation of r]{u)A'V{p{u) + V)-ncr on a specific boundary a between 
identified control volumes. They result from the variational formulation of the 
scheme and are viewed as fluxes between control volumes located at nodes and 
centers. 

Defining, for all k S A4 and u = {uh,Us)k.,s G Wxi, the diagonal matrix Mk(m) € 
by 

f , / ri{u„)+ri(u,') 

(Mk(ix)),,,, = V 2 itS-S, 

(0 otherwise, 

the systems (43) is equivalent to the following compact formulation: Vv G W-v, 

(45) f TTxiU^TT-Dvdx + Atn Y^ ^K.h{u'^) ■ 'Bf-{u'^)Sf-v = f 7rx)M"“^7rx)t’da;, 

Jq Jq 

where 

(46) B,^(w) := M„(m)A„M^(m), Vk G M, Vu G Wt>, 

is a symmetric semi-positive matrix since A,, and M^(m) have this property. 

2.3. Gradient flow interpretation for the scheme. The goal of this section 
is to transpose the formal variational structure pointed out in §1.2 to the discrete 
setting. A natural discretization of the manifold Wl consists in 


(47) 

leading to 

(48) 


3Jlx> = s u G Wt) 


TTxjudx = 


J Moda;| , 


TufXftx) = < ■« € Wxi 


TTxvdx = 0 


In order to define the discrete counterpart gx.u of the metric tensor g„ defined 
by (12)-(13), one needs a discrete counterpart of 

• the classical L^(D) scalar product: we will use 

{Wi,W2)l-^ / TTxWiTTxW2dx, VWi,W2GWx; 

Jn 
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the weighted H^{n) “scalar product” with weight ri{u): we use 

{WI,W2) ^ ^ S^Wi ■B,,{u)S,,W 2, yWi,W2 e W-D- 


kGM 


This allows to define the discrete metric tensor qti,u by: Vm S ITp, ywi,W 2 G 

TulH-p, 

(49) 

9v,uiwi,W2)= / tt-dWittt)4>2<^x = / 7r-D(/>i7rx)W2dai = ■ B^{u)6^(f)2, 


kGM 


where 4>i G Wx> solves the discrete counterpart of (13), that is 

(50) ^ d • B«,(M)d^'0 = / TT-DUijTT-Di/ida:, Vi/i G TTd. 


kGA4 


In this setting, we can define the semi-discrete in space gradient flow by 

(51) 9 v,u{dtu,w)= / TTx:{dtu)-Kxiv dx 

Jn 

= fl-D,'u(-V„£(M),i(j) = / TTvh{u) TT-pw da; = S^h{u) ■ B,„(m)(5^i;, 

where v solve the discrete elliptic problem 

^ d.i; • B K,{u)6i^xp = / TTxjui7r-p'0da;, \/ij} &W-d- 
kgot 

In order to recover (45) from (51), one applies the backward Euler scheme. 

Remark 2.2. In their seminal paper [65], Jordan, Kinderlehrer and Otto proposed 
to approximate the solution of gradient flows thanks to the minimizing movement 
scheme 

• C)(m, 16 ™“^) 


(52) 


u G argminl 


GSiti, ^ 2At 


+ 


<Sx)(m)|, 


where 0 denotes the distance on fOI-p induced by the metric tensor field gp. Several 
practical and theoretical difficulties arise when one aims at using (52). First of all, 
the Riemannian structure is formal, even in the continuous case. It is unclear if 
one can define rigorously a distance 5 if rjlff) — 0 even if g is concave (cf. [43, 81] j. 
But even if d is a distance, yielding a metric structure for dJl-E}, computing this 
distance is a complex problem we avoid by using an backward Euler scheme rather 
than (52). 

2.4. Main results. The first result we want to point out concerns the scheme for 
a fixed mesh. The following theorem states that the scheme (43) admits at least 
one solution, and justifies the free energy diminishing denomination for the scheme. 

Theorem 2.3. Let G IE™, then there exists (at least) one vector m” G 
solution to the system (43), and the following dissipation property holds: 


(53) 




kGM 


where Clp is defined by (36) and h(u^) = (f)K(RK), 1)3(11"))^ s defined by (38). 
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Since vP e W™ and since C W|,", Theorem 2.3 allows to define the iterated 
solution u — [u'^)i<n<N G to the scheme (43). 

The proof of Theorem 2.3 is contained in §3, together with some supplementary 
material that allows to carry out the convergence analysis when the discretization 
steps tend to 0. More precisely, we consider a sequence (T’m)m>i = 
of discretizations of 12 as introduced in §2.1.1, such that 

(54a) hj- = max hx —^ 0, 

T^Tm m-H-oo 

and such that there exists 0* > 0 and £* > 0 satisfying 
(54b) sup 9r^ < 9*, sup < £*, 

m>l m>l 

where 9j-^ and are defined by (21) and (24) respectively. 

Even though it can be avoided in some specific situations, we also do the following 
assumption, allowing to circumvent some technical difficulties: 

(54c) inf Cp„ = C* > 0. 

m>l 

This means that there is a minimum ratio of volume allocated to the cell centers 
and to the nodes in the mass lumping procedure. 

Concerning the time-discretizations, we consider a sequence (Atm)m>i 
cretizations of (0,tf) as prescribed in §2.1.2 : 

^t^n, = (Ati^^, . . . , At7V„,,m) ) Vm ^ 1. 

We assume that the time discretization step tends to 0, i.e., 

(54d) Atm = max At„ m —t 0. 

l<n<Nm ’ m—>-+00 

Theorem 2.4. Let (I?m, Atm)m be a sequence of discretizations of Qt^ satisfying 
Assumptions (54), and let At„) be a corresponding 

sequence of iterated discrete solutions, then 

7rp„,At„Mm —t u strongly in L^(Qti), 

m—^+oo 

where u is the unique weak solution to (1) in the sense of Definition 1. 

Proving the Theorem 2.4 is the purpose of §4. The practical implementation of 
the scheme (43) is discussed in §5, where we also give evidences of the efficiency of 
the scheme. 


3. Proof of Theorem 2.3 and additional estimates 

In order to ease the reading of the paper, several technical lemmas have been 
postponed to Appendix. 

3.1. One-step A priori estimates. 

Lemma 3.1. Let G bP™, and let tt" G he a solution to the scheme (43), 
then (53) holds. 
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Proof. Substituting v by h{u'^) = ( 18 ( 11 "))^ s defined by (38) in (45) yields 

(55) f — TT-pit"”^) 7r-p/i(it")da; + SK,h{u^) ■ Hf^h{u^) = 0. 

kGAI 

It follows from the convexity of F that 

r(a) — r(6) < (a — 6) (p(a) — p(l)), V a, 6 G M s.t. r(a), r(6) <+ 00 . 
Hence, using (44), one has 

f (tt-pm"' — TTpu"”^) 7rph.(M")da; 

Jn 


> 


= [ (tTpM" - TTpM” 

Jn 

[ (r(7rpM”) - r(7rpM 

Jn 


) (p(7rpM”) + 7rp(V)) da; 

^) + 7rp(M" — M"“^)7rpV) da; 

= €p(m") - £p(m"-i 


)• 


Using this inequality in (55) provides (53). □ 

Lemma 3.2. For all e > 0, there exists Cg G IR depending on e and p such that 

|u| < er(M)+ C',, VwGVFI,". 

Proof. Fix e > 0, then in view of Assumption (A2), the intermediate value theorem 
ensures the existence of Ue > 1 such that p{U() = p{l) + 1/e. Then for all u G Ip, 
one has 

pU pU 

T{u)= {p{a) - p{l))da = T{ue) + {p{a) - p{l))da. 

J 1 J Ue 

The function p being increasing, we deduce that 

r(M) > r(Ue) + {p{Ue) - p{l))\u -Ue\> r(Ue) + ^ (|u| - |Me|) , Vm G Ip. 

Lemma 3.2 follows with = |ite| — er(Me). □ 


Lemma 3.3. For all e > 0, there exists G M depending on e, rj and T such that 

ri{u) < er(u) + Ce, yu alp. 

Proof. The function u i-G tends to 0 as |it| —>■ 00 thanks to Assumption (A6). 
Let e > 0, then there exists > 0 such that 

|m| > Tj 0 < 77(11) < er(ii). 

The function 77 being continuous and nonnegative according to Assumption (Al), 
we know that 

0 < Ce := max 77 ( 11 ) < + 00 . 

u^[-rc,rc] 

The result of Lemma 3.3 follows. □ 

Lemma 3.4. There exist Ci and C 2 depending only on p, V and H such that 

(56) ^'Sp(it) + Ci< [ r{Tr-ou)dx < 2iBT>{u) + C 2 , Vu G FF™. 

^ Jn 

In particular, the discrete entropy functional £p is bounded from below uniformly 
w.r.t. the discretization P. 
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Proof. Recall that the discrete entropy functional C-p is defined by 

€p(m) = f (r(7rpM) + TrpVTTpit) da;, Vit G W|,". 

Jn 

Hence, one has 

(57) / r(7rpM)da: < €p(m) + ||7rpV||ioo(f2)||7rpM||ii(o), Vw G VF™. 

Jn 

Let e > 0 a parameter to be fixed later on. Thanks to Lemma 3.2, there exists a 
quantity depending only on p and e such that 

|w| < er(u) + Ce, Vu G Ip, 

ensuring that 

(58) ||7rpL6||Li(n) < e j r(7rpM)da; + C'emeas(f2), Vit G W|,”. 

Jn 

On the other hand. Assumption (A5) together with the definition (37) of V = 
(Vk, Vs)^ s ensure that 

lkpV|U.(0) < IIHIIoo. 

Setting e = 2 p 7 j| — in (58) and injecting the resulting estimate in (57) ends the 
proof of the second inequality of (56). The proof of the first inequality of (56) 
being similar, it is left to the reader. □ 

Lemma 3.5. There exists C depending only on A, fl, Op, fp, tp, rj, p and V such 
that, for all v = {vk.,Vs)k.,s G Wp^, one has 


H H H Ks'l p«.s(i^)(p(ii«)-p(^^s))^ 
kGA4sGV„ Vs'GV,, / 

<(7 11 + €p(t;) + ^ Sf,h{v) ■ B«,(i;)^k/i(«) I , 

V kGA 4 / 

where we have set ?7„_s(i;) = for all k € A4 and all s G Vk. 

Proof. Let v G C ILp", then it follows from the definition (38) of h{v) = 
(f)K(nK),f)s(ns))^^s € Wp that 

SM ■ + 2<5.V • Bk-«)^.V, Vk G M. 

It follows from Lemma A.2 that there exists C depending only on A, Op and ^p 
such that 

X! ( X! K.s'l) 7«.s('u) (p(^^k)-P( ns))^ < C^KP('y) ■ B«(i;)^kP(v), VkGM. 

sGV„ Vs'GV,. / 

Therefore, it only remains to prove that 

(60) ^ <5,V.Bk'*^)^«V< (Sp(i;) + (7, 

kGM 

for some C depending only on the prescribed data. Using Lemma A.3, we get that 


(61) S,^V ■B,^{v)S,^Y < Y max?7„,k^) Y ( 


kGM 


kGM 


(K - Us) . 


sGV., Vs'GV,, 
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It results from Lemma A.2 that for all k G AI, 


(62) E E Ks'l ^ 67 / VrV • AVrVdai < C|K;|A*||Vy||^. 

sGV„ Vs'GV^ / 

Denote by 77 ( 11 ) = {r]i^{v),f]^{v))^ ^ S Wxi the vector defined by 

rj^{v) = max r]{vk); max 77 ( 77 ') ) , Vsi'^) =0; yn G M, Vs S V, 

V s'6K / 


and remark that 


max?7K,s(L’) < 7 ?k(’w), Vk G M. 

SGVk 


Hence, we deduce using (62) in (61) that 

E • B,,(7;)5„V <C [ 7rMV{v)dx, 


kGM 


for some C depending only on Or, A, and ||VH||oo, the operator ttm being 
defined by (31). Let us now use Lemma A .8 to obtain that 


(63) 


E • B,(i;)^,V <C 7ri,77(i;)da7 


kGAI 


for some C depending only on the prescribed data, namely Or, A, fx), Qv and 
||VH||oo- Using Lemma 3.3, we know that for all e > 0, there exists depending 
only on e, 77 , L and meas(U) such that 


/ 7r-p77(i7)da; < e / r(7r-pi7)da; + Cg. 

Jn Jq 


Combining this result with Lemma 3.4 and (63), we deduce that for all e > 0, there 
exists Ce depending only on e, A, U, Or, C,v, V, P and V such that 

E • B,(7;)5,V < eC<Br{v) + Vv G 


kGM 

We obtain (60) by choosing e = i. This ends the proof of Lemma 3.5. 

o 


□ 


Lemma 3.6. Let ^ G W™, and let it" G he a solution to the scheme (43). 
There exist C\ and C 2 depending on At„, A, U, Or, (,- 0 , 0, P o,nd V such that 


E E ( E i<s'i)<.(p«)-p«))^ 

k&Ms&V^ \s'GV„ / 


<Ci(1 + (Sx,(m"-1)) <C 2 





Proof. Since m" is a solution of the scheme (43), the nonlinear discrete stability 
estimate (53) holds. Therefore, taking (53) into account in (59) yields 

E E ( E Kvl ) {p{0-p{u-)f < Cl (I + er{u--^)) 

KGMseV,, Vs'GV,, / 

for some Ci depending on the prescribed data. Then it only remains to use 
Lemma 3.4 to conclude the proof of Lemma 3.6. □ 
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3.2. Existence of a discrete solution. The scheme (43) can be rewritten in the 
form of a nonlinear system 

= 0r#ai+#v. 

In the case where p{0) = —oo, the function JF is continuous on , but not uni¬ 
formly continuous. The existence proof for a discrete solution we propose relies on 
a topological degree argument (see e.g. [73, 41]), whence we need to restrict the set 
of the possible rt” for recovering the uniform continuity by avoiding the singularity 
near 0. This is the purpose of the following lemma, which is an adaptation of [34, 
Lemma 3.10]. 

Lemma 3.7. Let G IL™ be such that 7rx)M"“^da; > 0 and let tt" be a 

solution to the scheme (43). Assume that p{0) = —oo, then there exists ev,At„ > 0 
depending on D, A, LI, rj, p, V, and £-p(m"'“^) such that 

K > e'D,At„, \/iy G MUV. 


Proof. First of all, remark that proving Lemma 3.7 is equivalent to proving that 
there exists Cxi.At„ > 0 such that 

(64) p{uZ)>-C-D,At„, MvGMUV. 

Because of the conservation of mass (44), we have 


n-pu'^dx = / TTpU^ ^da; > 0. 


Therefore, we can claim that there exists Vi G AiUV such that 

1 


(65) 


ul > 


TTxjit” ^da; > 0. 


meas(fl) Jq 

Let (z/f) G Ai U V he arbitrary, and {vq)q=o,...,i be a path from Vi to Vf, i.e. 

• vq = Vi, VI = V{, and t-p ^ if p ^ g; 

• for all g e {0,..., — 1}, one has: 


Vq G M Vq+l G and Vq G V Vq+i G 

Let q G {0,... — 1} 

It follows from Lemma 3.6 that there exists Cp^At„ depending on V, At„ A, LI, 
ri, p, V, and C:-d(m"“^) such that 

K,eM sGV„ 

This ensures in particular that 
e-i 

(66) < Cv,At„, 

q=0 

where we have set = Vk,s if Wq^i^q+i} = {k,s}. 

We can now prove (64) thanks to an induction along the path. Assume that 

Kg > ex>.At„ for some e-D,Az„ > 0, whence - ^T’.At„ > 0- Tfien it 

follows from (66) that 

^ ~^T>,At„ '*^",+1 — ^x>.Az„ > 0- 


C- 


V,At„ 


'^V.At„ 


P{Kgj>P{Kg) 
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We conclude as in [34, Lemma 3.10] thanks to the finite number of possible paths. 

□ 


Thanks to Lemma 3.7, one can apply the same strategy as in [34] for proving 
the existence of a solution to the scheme (43). 

Proposition 3.8. Let € W™, then there exists (at least) one vector tt" G 

solution to the system (43). 

Proof. As in [34, Proposition 3.11], the proof relies on a topological degree argument 
(cf. [73, 41]) applied twice. More precisely, start from the parametrized nonlinear 
problem that consists in looking for solution to 


(67) 


TTX) 




ttvvAx + 7 XI XI - f)s«’^)) {vk - v'f) 

sGV,, s'GV„ 


+ (1 - 7 ) X ( X - fls«’^)) {v^ - Ts) = 0, Vn e W-D- 

sGV„ Vs'GV,, / 


For 7 = 0, the corresponding scheme is monotone, hence the nonlinear system (67) 
admits a unique solution 14 ”’° G and its corresponding topological degree is 

equal to 1 (see for instance [49] for the application of this argument to the case of 
a pure hyperbolic equation). In the case where lim„N^o7'('^) = — 00 , one can prove 
as in Lemma 3.7 that any solution m"’"'' G to (67) satisfies 

u}'^>evAtr., V/3eMUV 


for some e-v.Atr^ > 0 not depending on 7 . The convex subset on which one looks for 
a solution can be restricted to the subset defined by 


K := \^uG 


U 0 > and €-d(m) < C:-d(m” ^) + l|, 


the corresponding topological degree being still equal to 1. Note that the bound 
up > must be removed if p(0) is finite. This ensures the existence of at least 

one solution to the nonlinear system (67) when 7 is equal to 1. 

Starting from the system (67) with 7 = 1 , one defines a second homotopy 
parametrized by /r S [0,1] to get (45). More precisely (the superscript 7 = 1 
has been removed for clarity), we set 


ri^ : ui-G- 1 + fi{ri{u) - 1), G [0,1], 

so that ?7° is constant equal to 1 and 77^ = 77. Define G as a solution to 


_ 7i"-l ' 

(68) / TTX) ( -- ) TTvvdx 

Jn 


Atr 


+ X X - f)s«’'')) (^'« - vi) = 0, 

sGV„ s'GV„ 

where v G VFp is arbitrary, and where 

n.M ^ 77^(0 + 77^(0 

''K.S' 


2 


Vk G M, Vs G Vk. 










22 


CLEMENT GANCES AND CINDY GUICHARD 


A priori estimates similar to those derived previously in the paper ensure that 
the solutions to ( 68 ) cannot belong to dK whatever the value of /r S [0,1]. The 
existence of a solution to the scheme (43) follows. □ 

3.3. Multistep a priori estimates. As a byproduct of the existence of a discrete 
solution It” for all n G {1,..., N}, we can now derive a priori estimates on functions 
reconstructed thanks to the discrete solution u G Wx),At- 

The first estimate we get is obtained by summing Ineq. (53) w.r.t. n, and by 
using the positivity of the dissipation. This provides 

(69) max £x)('i‘"') < <S-d('W°) < S(wo) + C <C, 

nG{l,...,IV} 

where C only depends on V, uq, p and 12 thanks to Lemma 2.1. Since the discrete 
entropy functional (£x> is bounded from below by a quantity depending only on p, 
V and n (cf. Lemma 3.4), we deduce also from the summation of (53) w.r.t. n 
that there exists C depending only on p, V, 12, and uq (but not on V) such that 

N 

(70) ^ ^ d^u^) ■ < C. 

n—1 kGAI 

Mimicking the proof of Lemma 3.5, this yields 

N 

(71) ^ ^ 5,p{u^) ■ B.(m”)5«p(m") < C 

n—1 kGA4 

for some quantity C depending on A, 12, 6 * 7 -, C,v, ly, P, V and tf. 

The following lemma is a direct consequence of Estimate (69) and Lemma 3.4. 
Its detailed proof is left to the reader. 

Lemma 3.9. There exists C depending only on 12, p, V, uq and 12 (but not on the 
discretization) such that 

|jr(7rx)^AtM)|l7^oo((o,tf);Li(n)) ^ C*. 


The introduction of the Kirchhoff transform was avoided in our scheme. Its 
extension to complex problems (like e.g., systems, problems with hysteresis) in 
therefore easier. However, the (semi-)Kirchhoff transform ( defined by (7) is useful 
for carrying the analysis out. The purpose of the following lemma is to provide a 
discrete L^((0, T); 77^(12)) estimate on (,{u). 


Lemma 3.10. There exists C > Odepending only on A, p, V, 9j- and Id, 12, 2f, 
uo such that 

AVr,At^(M) • Vr,AtC(M)da; < C. 

Qtf 

Proof. Since the function rj was assumed to be nondecreasing, see Assumption (Al). 
We know that for all interval [a,b] C Ip, maXf,g[a,6] 17 (c) = max{r 7 (a), 77 ( 6 )}, hence, 
denoting by X"g the interval with extremities u” and u", we obtain that 


(72) 


i7ks ^ A i^iax 77 (c), Vk G M, Vs G Vk. 


- 2 ciij 3 

The definition (7) of the function ( implies that 

(?«) - m)f < (max 77(c)) {p{u:)-p[ul)y 
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whence we obtain that for all k & M., 

\M,iu-)SMnnf = E {PK)-P{u:)f 

s6V^ 


Using that 

V ■ A^v > w ■ AkW, s.t. |ii|^ > Cond 2 (AK)|ip|^, 

and that = M„(m)A„Mk(m), we get that for all k G A4, 




1 


2Cond2(A„ 

1 

2Cond2(A„ 


^ [ AVr?(M’") • VrC(M”)da;. 


Thanks to Lemma A.l stated in appendix, we know that C > 0 depending only 
on A,0j- and £xi such that Cond 2 (AK) < C, for all k G M, so that: Vk G M, 
Vn e A}, 


(73) 


Vre(M") • AVre(w")clai < Cd^u^) ■ 


In order to conclude the proof, it only remains to multiply (73) by At„, to sum 
over K G M. and n G {1,..., A}, and finally to use (71). □ 


Combining Estimate (71) and Lemma A.2 yields the following lemma, whose 
complete proof is left to the reader. 

Lemma 3.11. There exists C depending only on A, 11, Oq-, (x>, P> P, V o,nd if 
such that 

E E ( 

n—l kGAI sGVk X! 



E 

'^Vk 


a:j]v:jp{u:)-p{u:)f<c. 


4. Proof of Theorem 2.4 

In what follows, we consider a sequence (Dm, Atm)m>i of discretizations of Qtf 
such that (54) holds. In order to prove the convergence of the reconstructed discrete 
solution towards the weak solution of (1) as m tends to oo, we adopt 

the classical strategy that consists in showing first that the family (7rx)„,,At„'Um)„j>]^ 
is precompact in L^{Qtf) (this is the purpose of §4.1), then to identify in §4.2 the 
limit as a weak solution of (1) in the sense of Definition 1. 

As a direct consequence of Theorem 2.3, one knows that the scheme admits 
a solution Um = rm^s m) that, thanks to the regularity assumptions (54b)- 
(54c) on the discretization and thanks to Lemmas 3.9, 3.10 and 3.11, satisfies the 
following uniform estimates w.r.t. m: 

Atmr(Mm)|lioo((o,tf);Li(n)) — 


Vr„.At„C(Mm) • AVr„.At,„C(“m)da;dt < C, 
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/ 

(76) ^ EK«'I 

71—1 SGVk \sGVk 

where C may depend on the data of the continuous problem, and on the discretiza- 
tion regularity factors 0 *, i* and C* but not on m. 

4.1. Compactness properties of the discrete solutions. 

Lemma 4.1. Let {Dm, ^tm) be a sequence of discretizations of Qt^ satisfying As¬ 
sumptions (54), there exists C depending only on A, 9*, £*, f2, t{, p and uq such 
that, for all m> 1, one has 

IKTr„,At,„^('Wm,)|| 2 , 2 (('Q < C and ||TX)„,At„'f('Um)||^ 2 ('Qj^) < C. 

Proof. It follows from the Estimate (75) and from Lemma A.5 that for all m > 1, 

(77) ~ '^Tm.At„^('l*m) II j;,2((o,tf);Li(n)) 

< meas(n)^/^ ||7rx)„,At„'?(Mm) - 7rr,^,At„C(t*m)llL2(Q^^) < C 

for some C depending only on A, Ll, t{, £*, 9*, p, and uq (but not on m). 
Moreover, it follows from Assumption (8) that 

||7rx),„,At„$(U'm)||ioo((o ij).ii(Q)) < c* < C 

thanks to the estimate (74). Combining this inequality with (77) provides that 

IKTT„,Atm^('^m)llLi(Q,j) ^ C 

whence the sequence (7rr„,At„C(“m))„>i is bounded in L‘^{{0,t{)-,W^’^{fl)) thanks 
to (75). A classical bootstrap argument using Sobolev inequalities allows to claim 
that it is bounded in L^((0, tf); H^{Ll)), thus in particular in L?{Qtf). One concludes 
that (7r7-^_At„,C(Mm))„j>i is also bounded in L'^{Qtf) thanks to (77). □ 

Remark 4.2. Art alternative way to prove the key-point of Lemma 4-1, namely 

||7rr„.At„^(Wm)||j;,2((o < C, 

would consist in using [63, Lemma A.l], that shows that 

V G L\il) andVfiv) G L'^inf ^{v) G H\n). 

As a consequence of Lemma 4.1, we know that the sequence (7r7-^.At,„C('*^m))j,j>]^ 
is relatively compact for the L^((0, if); 77^(0))-weak topology. Moreover, the space 
being locally compact in L^(0), a uniform information on the time translates 
of 7r7-^_At„C(wm) will provide the relative compactness of (7r7-^_At„C('Um))„^>]^ in 
the L^((5tj)-strong topology (see e.g. [91]). Such a uniform time-translate estimate 
can be obtained by using directly the numerical scheme (see e.g. [51, 34]). One 
can also make use of black-boxes like e.g. [6, 8]. Note that the result of [60] does 
not apply here because of the degeneracy of the problem. We do not provide the 
proof of next proposition here, since a suitable black-box will be contained in the 
forthcoming contribution [8]. A more classical but calculative possibility would 
consist in mimicking the proof of [34, Lemmas 4.3 and 4.5]. 


v:Ap«J-pKm)Y<c, 
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Proposition 4.3. Let (Vm, ^tm) be a sequence of discretizations of Qt^ satisfying 
Assumptions (54), and let be the corresponding sequence of solutions to 

the scheme (43). Then there exists a measurable function u : Qt^ —>■ K with f{u) € 
L^((0,t{); H^(fl)) such that, up to a subsequence, one has 

(78) U a.e. inQt^. 

m—¥oo 


Corollary 4.4. Keeping the assumption and notations of Proposition 4-3, one has 

u strongly in L^{Qtf). 


Proof. As a result of Proposition 4.3, the almost everywhere convergence prop¬ 
erty (78) holds. On the other hand, it follows from Assumption (A2), more pre¬ 
cisely from the fact that lim„_>oo p{u) = -l-oo that the function P defined by (4) is 
superlinear, i.e.. 


lim 


r(w) 


= - 1 - 00 . 


«—>- 1-00 u 

Therefore, Estimate (74) implies that (7i'x),„,At„t*m)m>i is uniformly equi-integrable. 
Hence we can apply Vitali’s convergence theorem to conclude the proof of Corol¬ 
lary 4.4. □ 


Lemma 4.5. Under the assumptions of Proposition 4-3, one has 

7rr„.At„C(Mm) —!■ C(w) weakly in L^{{Q,ti)-,H^{Ui)), 

m—>-oo 

where u is the solution exhibited in Proposition 4-3. 

Proof. Thanks to Lemma 4.1, the sequence (V7-^^At„'f ('Urri))m>i i® uniformly bounded 
in Therefore, there exists S € L^((0,tf);i7^(H)) such that 

7rr„,At„C(«m) —!• S weakly in L^((0,tf); 77^(H)). 

m—^oo 

But in view of Proposition 4.3 and of the continuity of we know that 
^(7r-D,„,At„Wm) = 7rx)„,At,„'?(Mm) — f{u) a.e. in Qtj. 

m—>oo 

Since and 7rx)„,,At„C('*^m) have the same limit (cf. Lemma A. 5), we 

get that E = f{u). □ 

Lemma 4.6. Let u he the limit value of obtained in Proposition 4-3, 

then 

(79) 7r-D,„,At„?y(Mm) — viu) strongly in L^{Qt^), 

771—>-00 

and 

(80) viu) strongly in L^{Qtf). 

771—>-00 

Proof. Let us first establish (79). Thanks to the entropy estimate (74), we know 
that the sequence (7r-p„,At„r(Mm)),„>;^ is uniformly bounded in L°°((0, tf); L^(H)), 
thus in L^(Qtj). Then Assumption (9) allows to use the de la Vallee-Poussin theo¬ 
rem to claim that is uniformly equi-integrable on Qt^. More¬ 

over, the continuity of ry and Proposition 4.3 provide that 

(7rx),„,At„iy(Mm))„,>i —!• viu) a.e. inQtj. 

— 771—>-00 

Therefore, we obtain (79) by applying Vitali’s theorem. 
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Let US now prove (80) by proving that 7rx)„,At„i?(wm) and 7rx„,At„?7(itm) (that 
is uniformly equi-integrable for the same reasons as is) have the 

same limit r]{u) as m tends to oo. It follows from a combination of (75) with 
Lemma A.5 that, still up to an unlabeled subsequence, 

(81) 7’'X)„,At„'f(Mm) - 7rvvt„,At^'f(Mm) —^ 0 a.e. inQjj. 

m—>-oo 


Since the function o is assumed to be uniformly continuous (cf. (10)), it 
admits a non-decreasing modulus of continuity w € M+) with ti7(0) = 0 such 

that, for all v,v in the range of 


(82) 


y/r]o^ i(u) - y/r]o ^ i (u) 


< w 


u-u|). 


m) 7rA4„,At^C(Wm)|)- 


SO that 

-\/7r-D„,At„?7(M m) m) 

Therefore, it follows from (81) that 

\/7rx)„,At„i?(Mm) - \/7rAi„,,At„??(Mm) —0 a.e. inQtf. 

V V m—>-oo 


Thus (7r-D„,At„i7(“m))^>i and ( 717 ^ 4 ^,At„??(Mm))„>i share the same limit. □ 


4.2. Identification of the limit as a weak solution. 


Proposition 4.7. Letu be a limit value of the sequence (7r-p^^At„'*.tTn)„>i exhibited 
in Proposition 4-3, then u is the unique weak solution to the problem (1) in the sense 
of Definition 1. 


Proof. In order to check that u is a weak solution, it only remains to check that the 
weak formulation (11) holds. Let ip G Cp°{D. x [0,tf)), then, for all m > 1, for all 
p G A4mU Vm and all n G {0, ..., N^}, we denote by Atm = (Ati,™, • ■ •, AtN^,m), 

by tn.m = Er=l by Ip^ = 1p{xis,tn,m), by Xpm = (tpp) G W-D^, 

and by ipm = {4’m)o<n<N^ ^ lAD„,Ai^- Note that since ip{-,tf) = 0, one has 
ipm^ = 0 for all m > 1. 

Setting V = xp^^ in (45) and summing over n G {1,. ■ ■, Nm} leads after a 
classical reorganization of the sums [51] to 

(83) Am + Bm + Cm + Dm = 0, 


where we have set 


Am — ^ ^ ^tyi m I '^T>jnXL.yy^{x')lTx>jn [ 
n=l 3 ^ V 

Bm=- '^v^u^(x)Trx,^-ip°(x)dx, 
Jn 


^^n.m 


{x)dx^ 


Nrr 


Cm=Y.At^,m Y. ■ B.(n” )d«V’[r\ 


n—1 

Nm 


K^Mr. 


Dm=Y^t^,rn Y ^ m ' B ,{u^)S ,xPPPf\ 


n—1 


KGMr. 


and Vm — (14(3;^), P(a;s))„g 7 V(„,sGV„ ■ 
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The regularity of "0 yields 


Nm 

n—l 


^"-1 


-dt'ip uniformly on Qtj 


where tn,m = Er=i so that, using Corollary 4.4, one gets 


(84) 


lim Am = - 

m—>-oo 


udtip dxdt. 


Qt, 


The function TT-p^Umix) tends strongly in L^(Sl) towards uq and 7r-p„'i/’° con¬ 
verges uniformly towards ip{-, 0) as m tends to -l-oo, leading to 

(85) lim Bm = — / uo(a;)^(x, 0) da;. 

m-)-oo 

We split the term Cm into three parts 

(86) Cm = Ci^m + C2,m + C^^m, TO > 1, 

where, setting = xp^m^ '^m = for all n G {!,...,iV^}, and 


= (jpl 


0 <n<Nrr. 


€ one has 


C^l,m — II 7^ ,Ai tjt,^, Airn 


N„ 


C2,m=^^tn,m T T V^s (piO - PK)) <, 3 ' 

n—l K^Aiyn sGVk s'GVk 

X - VviK)) ^ 

Nrr, 

C3,m='^Atn,m T iVWs{p {0 “ P«)) “ (^K) “ ?«))) 

X E 


n—l 


KeM„ 


sev^ 


s'ev^ 

Thanks to Lemma 4.6, we know that 

7r7w„,At„,\/??(Mm) —vTU Strongly in U((3tf). 

m—^oo 

Hence, it follows from the weak convergence in L?{Qtf) of Vr,^,At„C(wm) towards 
V^(u) (cf. Lemma 4.5) and from the uniform convergence of ^towards 
Vi/> as TO tends to -l-oo (see for instance [40, Theorem 16.1]) that 


(87) 


lim Cim= // \/77(u)V^(u) • AV'0 da;dt. 

’ JJot, 


Let us focus now of (72,™• Using the inequality ah < ea^ + one gets that 


(72,™ < eC -I- 


r^n 
0^2,m 

YU’ 


( 88 ) 


Ve > 0, 
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where 


Nrr 


E E E Ks'l 


n—l 

Nrr. 


KG:Mm. SGVk Xs'^Vk 


n=l kGVW™ sGV„ Vs'GV^ / 

We deduce from Estimate (76) that 

(89) < C, Vm > 1, 

for some C depending only on uq, p, ft, A, 9* and i*. 
Define ^ by 

V” = 0, VsGV„, 


(90) 


fJ-K = maxsev„ 


Vv{us) - VviK) 


Vk G Mm, 


Vn G {0,.. .,Nm}- 


The definition (42) of ? 7 " g implies that 

VviK) 


^ < /x”, Vk G TWm, Vs G V„. 

Therefore, we get that 

C'l^<^At„,m E E f E K^'I) 

n=l kGXV(„ sGV,, Vs'GV,^ / 

Then thanks to Lemma A.2, there exists C depending only on A, 6* and £* such 
that 


C" < C 


Qt. 


(7’‘A1„,At,„Mm) ’ AV 7 -,„,At,„V’md£cdt, Vm > 1. 


Since 7r2n^_A.t„Mm converges to 0 strongly in L'^{Qti) as m tends to oo (this is 
the purpose of Lemma 4.8 hereafter), and since remains bounded in 

L°°{Qtt) uniformly w.r.t. m > 1, one gets that 

(91) lim C'l^ = 0. 

m—^oo 

Therefore, it follows from (88)-(91) that limsup^_,.go C^^m < Ce for arbitrary small 
values of e > 0, whence 

(92) lim C 2 .™ = 0. 

m—^oo 

As a preliminary before considering let us set, for all k G Mm, all s G V„, 

and all n G {1,..., Nm}, 

2 


' ao-ao 

^.s = i VMw”) -p(<) 
MO 


if w" = 


Thanks to the mean value theorem, we can claim that, for all k G Aim, all s G Vk 
and all n G Nm}, there exists w" g G = [min('u", m"), max(u", u”)] such 

that 77 k s = vMk s)- la particular, this ensures that 


<M, 


Vk G Mm, Vs G Vk 
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where was defined by (90). Using moreover that < 277 ^^, one gets that 


Nrr 


C3,m < 2 ^ E E V<sipi<) - pK)) e <s' (c-' - V'r'). 

71—1 sGVk s'eV„, 


Cauchy-Schwarz inequality and Estimate (76) yield 

G,™ < c f E E i E K^'i) -C"')' 

KkGM^ sGV^ \s'GV„ / 

Lemma A.2, and the regularity of tp provide that 

C'a.m ^ C \\T^Mm-^trrit^m\\L^{Qtf) ■ 

Applying Lemma 4.8, we get 

(93) lim Ca,™ = 0. 

m—xx) 

Putting (86) together with (87), (92), and (93), one gets that 

(94) lim C^ = II ■ AV^ da;dt. 

JjQt^ 


1/2 


Now, we focus on the term Dm that can be decomposed into 
(95) 

where we have set 


— 7)i m 4“ -^2,m 4“ -^3,mi V?TZ ^ 1, 


7^1,m — 


7r;w™,At„f?(Mm)Vr,^V„ • AVr^.At^V’mda^dt, 


Qt 

Nrr 


7?2,m — 2 'y ^^n,m ^ ^ ^ ^ ' ®s,s' (Lt 14) 

n—1 K^Aim sGVk 

X {\f^' + VviK)) - i’sA), 

Nm 

Dz,m =2 y ^tn,m ^ ^ ^ ^ ^ ®s,s' (•\/''7 k,s + a/ 1?(m")) (Vk ~ 14) 

n—1 KieAtm sev^ 


It follows from Lemma 4.6, from the uniform convergence of Vt-^V^ towards 
VE and of towards Vi/j as m tends to +oo that 


(96) 


lim Dim= / 77('u)VF • AV'0da;dt. 

’ Jo,. 


Let e > 0, using again the inequality |a6| < ear + we obtain that 


|A’2,m| < 


(97) 


Vm > 1, 
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where we have set 

Nrr 


D'2,m=^^in,m XI XI ( H 1 (\/^+\A(^) (C ^)^ 

n —1 sev^ Vs^ev^ / 

XAt„,„i X X fx Ks'l^ (y^-vX^) (i4-v;)L 


D — 

^2,m — 

n—1 KGA4m sev^ \s^eV^ 

Define rj^ = e Wx.„,At„ by 


Vl = max 


then one has 


whence 

^ 


(? 7 «),max<), ?7” = 0> Vk e Vs G V„, 

\ sGV„ / 

/? 7 ^+ VviK)) <47?", Vk G Vs G V^, 


Qt 


7r;w™.At™'^m AVr„,At„'0m • Vr„,At„'*/’mda;dt, 


Vm > 1 


thanks to Lemma A.2. Using Lemma A.8, we know that there exists C depending 
on the data of the continuous problem and of the regularity factors 6*, £* and C,* 
such that 

hMrr^^t^flJl Ll(Qj^) < > 1, 

while the regularity of tp ensures that 


|Vr„,At„'0,, 


<C, 


Vm > 1. 


Therefore, there exists C depending only on the data of the continuous problem 
and the regularity factors 9*, £* and (* such that 

(98) <C, Vm > 1. 

The term D 2 ^ can be studied as was, leading to 

(99) hm D'l^ = 0, 
whence, taking (98)-(99) into account in (97), one gets that 


( 100 ) 


lim D 2 ,Tn = 0. 


Reproducing the calculations carried out for dealing with Z? 2 ,m allows to show that 
(101) lim Ds rn = 0. 

m—^oo 

Combining (95)-(96) and (lOO)-(lOl), we obtain that 


( 102 ) 


lim Drr, = 

m—>oo 


/ / r]{u)A'W ■ Vpjdxdt. 
■I-IQh 


Finally, it follows from (83), (84)-(85), (94) and (102) that the limit u of the 
discrete reconstructions (7rx)„,,Atm'*^m)m>i solution to the problem (1) in 

the sense of Definition 1 . □ 


Lemma 4.8. Let = (m",m") 


K ’ JK,S,n 


G Wi 




be defined by (90), then 




0 strongly in 
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Proof. Using (a — b)^ < 2 (a^ + 6 ^), one gets that 


< 2 


??«) +max? 7 «) < 2 

SGVk, 




s6V„ 




for all K S and all n G {1,..., Nm}- Using {54b)-(54c), which ensure that 

c* c* 

rriK. > —meas(K) and rus > -^meas(K), 

a di* 

we deduce that there exists C depending on d, i* and such that 

< CTTTi^,A.trnV{Um), Vto > 1. 

As a particular consequence of Lemma 4.6, we know that (7rx),„,At„i7(itm)),„>i is 
uniformly equi-integrable, whence 

(103) uniformly L^-equi-integrable. 

Let us introduce Wm = (w",s,n G defined for all k G Mm, all 

s G Vm and all n G {1,..., Nm} by 


(104) Ws” = 0, w” =max|^(M”)-C«)|. 

sGV„ 

It follows from a straightforward generalization of Lemma A.9 and from esti¬ 
mate (75) that converges strongly in Lf{Qtf) towards 0. Therefore, 

up to an unlabeled subsequence, it converges almost everywhere. As a consequence, 


(105) 


'^M,r..^trr,4’{Wm) - S’ 0 U.e. lu Qjj 

771—^00 


for all continuous function (f : M+ —>■ M such that ^( 0 ) = 0 . 

It follows from the definition of /.t^ that 

H]}<w{w2), Vk G Aim, Vn G {!,..., Am}, 
where zn is a modulus of continuity of -^770 (cf. (82)). Hence, we obtain that 

0 ^ , Ai„, ^(^Um ) ■ 

Thanks to (105), we obtain that 


(106) TTAi^.At^Mm —^ 0 a.e. inQtj. 

m—>-oo 

In order to conclude, it only remains to remark that (103) and (106) allow us to 
use Vitali’s convergence theorem. □ 


5. Numerical implementation and results 

This section is devoted to the numerical resolution of the nonlinear system (43). 
First, we discuss in §5.1 the strategy that we used for solving the nonlinear sys¬ 
tem (43). Then we present in §5.2 two 2-dimensional cases with analytical solutions 
in order to illustrate the numerical convergence of the method. 
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5.1. Newton method, Schur complement and time-step adaptation. The 

nonlinear system (43) obtained at each time step is solved by a Newton-Raphson al¬ 
gorithm. Given S Wp, this leads to the computation of a sequence (it"’*) C 
Wti such that it" = limi_,.oo it"’® is a solution to (43). The variation of the discrete 
unknowns between two Newton-Raphson algorithm iterations is denoted as follows, 


dit"’® = ( d<’®, d<’® ) 


sev,KeM 


= it"’®+i - It"’®, 


Vi > 0. 


Let us briefly detail the practical implementation of the iterative procedure allowing 
to deduce it" from it"“^. 


(1) In the case where p(0) is finite, the initial guess for the Newton algorithm 
is, as usual, taken as it"’° = (u"”^, it"“^)^ ^ for all s € V, tt € A4. In the 
singular case p(0) = —oo, it was proved in Lemma 3.7 that the solution 
It" = (m",u")k,s of (43) is such that min^g^viuv > 0. Therefore, we can 
initialize the Newton algorithm by 

it"’° = (max(e,u""^) ,max(e,M""^))^_^. 


In the computations, we fixed e = 10“^°. 

(2) The Newton-Raphson algorithm iterations are done until a convergence 
criterion on the L°°(f 2 ) norm of the variation of the discrete unknowns is 
reached or until the maximum number of iterations is reached. At each 
iteration, the Jacobian matrix resulting of (43) is computed and has the 
following block structure 




dit 


n,i _ 



where the sub-matrices have the following sizes: A G 0 B G 

R#^ 0 R#^, C G and D e The sub-vectors at 

the right hand side have thus the following sizes: bi G and b 2 G 
The dependence of the sub-matrices and the sub-vectors w.r.t. n and i was 
not highlighted here for the ease of notations. A main characteristic of this 
block structure is that the block D is a non singular diagonal matrix, thus 
the Schur complement can be easily computed without fill-in to eliminate 
the variation of the cell unknowns. This allow to reduce the linear system 
to the variation of the vertices unknowns as is usual when using the VAG 
scheme. The resulting linear system that we have to solve in order to obtain 
the variation of the vertices unknowns is given by, 

(107) (A - BD-iC)(d<’®),ev = bi - BD’^ba, 


and then the variation of the cell unknowns can be easily deduced by the 
matrix-vector product below. 


{dur).GM = D-i(b2 - C(d<’®),ev). 


As for the initial step, we have to take into account the singular case at 
each Newton-Raphson iteration by, 

m"’®+i = max(it"’® -k dit"’® , e). 

(3) If the Newton-Raphson algorithm stops before the maximum number of 
iterations is reached, the next time iteration is proceeded by increasing the 
time step. Otherwise, the current time iteration is recomputed by reducing 
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the time step. The time step is bounded by a maximum value denoted 
^^max- A maximum number of convergence failures of the nonlinear meth¬ 
ods is imposed in order to abort the simulation in case of a non-convergence. 

5.2. Definitions of the test-cases and numerical results. We present here 
four 2-dimensional numerical cases where is the unit square. The space domain 
is discretized by using meshes obtained from a benchmark on anisotropic diffusion 
problem [61]. In the following numerical experiments, the tensor is defined by 



where lx and ly are chosen constant in fl, and the exterior potential is defined by 
V{x) = —g • X for all a; e n where g = ( 5 , 0 )* with g G K_|_. The weights of 
the VAG scheme defined in (27) are defined by for all k G At, s G V^. 

We refer to [54, 25] for a discussion on the mass distribution for heterogeneous 
problems. The linear solver applied to solve (107) is a home-made direct solver 
using a gaussian elimination with an optimal reordering. 

In some of the test cases presented hereafter, Dirichlet boundary conditions 
are considered instead of no-flux boundary conditions. This allows to construct 
analytical solutions to the continuous problem and to perform a convergence study. 
Even though it has not been done in this paper, the convergence proof for the scheme 
can still be carried out when (sufficiently regular) Dirichlet boundary conditions are 
considered. However, the gradient flow structure is destroyed and the free energy 
might not be decreasing anymore in this case. 

Errors are computed in the classical discrete L‘^{Qti), and 

norms. All the results are presented in the Tables below. Each table provides the 
mesh size h, the initial and maximum time steps, the discrete errors, their associated 
convergence rate and the minimum value of the discrete solution. It also contains 
the total (integrated over time) number of Newton-Raphson iterations needed to 
compute the solution as a indicator of the cost of the numerical method. 

5.2.1. Test 1: Linear Fokker-Planck equation with no-flux boundary condition. This 
first test case matches with the problem defined by ( 1 ) with the functions giu) = u 
on M_|_ and p{u) = log(u), and with the gravitational potential V(x, y) = —gx where 
the constant g is fixed to 1. Setting g = ( 5 ,0)^, the problem (1) leads to the linear 


equation 


(108) 


- V-(A (Vu - ug) ) = 0 in Qtf 


We compare the results obtained with the nonlinear scheme (43) with those ob¬ 
tained using the definition of the fluxes 


(109) 



s'GV„ 


s'GV„ 


Vk G A4, Vs G V„ 


instead of (43c). The resulting scheme is called the linear scheme. The numerical 
convergence of both schemes has been compared on the following analytical solution 
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(built from a 1-dimensional case): 


( 110 ) 


u{x, y, t) = exp 



^7rcos(7ra:) 



-I- Trexp 




y{{x,y),t) €fix (0,t{), 


where a = lx ( 71 ^ + ^)- This function satisfies the homogeneous Neumann boundary 
condition and the property u{x,y,t) > 0 for all {x,y,t) G Qt^. 

In order to make a numerical convergence study, we have used a family of tri¬ 
angular meshes. These triangle meshes show no symmetry which could artificially 
increase the convergence rate. This family of meshes is built through the same pat¬ 
tern, which is reproduced at different scales: the first (coarsest) mesh and the third 
mesh are shown by Figure 4. Although the analytical solution is one-dimensional 
and the permeability tensor is diagonal, the discrete problem is really 2D because 
of the non-structured grids. The 2D aspect of the problem is amplified by the 
choice of a stronger diffusion in the transversal direction. For the tests on triangu- 



Figure 4. First and third mesh used in the numerical examples. 


lar grids, the final time t{ has been chosen to 0.25 and an anisotropic tensor has 
been consider: lx = 1 and ly = 10. 


h 

#v 

^^init 

^^max 

erri,2 

rate 

err/,1 

rate 


rate 

Ujnin 

:?^Newton 

0.250 

37 

0.001 

0.01024 

0.196E-01 


0.754E-02 

- 

0 . 216 E +00 


0.022 

204 

0.125 

129 

0.00025 

0.00256 

0.512E-02 

1.935 

0.178E-02 

2.084 

0 . 600 E -01 

1.848 

0.004 

456 

0.063 

481 

0.00006 

0.00064 

0.129E-02 

1.986 

0.430E-03 

2.050 

0.157E-01 

1.931 

0.001 

1307 

0.031 

1857 

0.00002 

0.00016 

0.324E-03 

1.997 

0.107E-03 

2.007 

0.473E-02 

1.734 

0.000 

3935 


Table 1. Triangles. Nonlinear scheme (43). 


h 

#v 

Aiinit 

Atmax 

err^2 

rate 

err^^i 

rate 


rate 

Umin 

:?^Newton 

0.250 

37 

0.001 

0.01024 

0.187E-01 


0.708E-02 


0.225E+00 


-0.155 

33 

0.125 

129 

0.00025 

0.00256 

0.469E-02 

1.993 

0.165E-02 

2.100 

0.786E-01 

1.515 

-0.046 

106 

0.063 

481 

0.00006 

0.00064 

0.117E-02 

1.999 

0.406E-03 

2.023 

0.228E-01 

1.784 

-0.012 

400 

0.031 

1857 

0.00002 

0.00016 

0.293E-03 

1.999 

0.102E-03 

1.999 

0.611E-02 

1.901 

-0.003 

1570 


Table 2. Triangles. Linear scheme, fluxes defined by (109). 


Let us first observe that the numerical order of convergence is close to 2 for both 
schemes. The nonlinear scheme is of course more expensive than the linear one but 
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it preserves the positivity of the solution, unlike the linear scheme. This numerical 
behavior is a verification of the theoretical result mentioned in the Lemma 3.7. In 
the linear case, the number of Newton-Raphson iterations is equal to the number 
of time steps. On the finest mesh, the ratio of the number of Newton iterations 
between the nonlinear and the linear schemes is about 2.5. It seems to be acceptable 
in cases where preserving the positivity is mandatory. 

Now, in order to exhibit the ability of the VAG scheme to deal with general 
meshes, the same test case has been applied on a so-called Kershaw grid (cf. Figure 
5). Instead of an irrelevant numerical convergence study — it is difficult to define 
a refinement factor for this type of grids —, we aim to give an evidence that the 
scheme is free energy diminishing (thus positivity preserving) and that the long¬ 
time behavior of the continuous problem is preserved at the discrete level by the 
scheme. The final time tf has been chosen to 250 and an anisotropic tensor has 
been consider: = 0.001 and ly = 1. The results are listed on the Table 3 and 

we can check again that the nonlinear scheme is positivity preserving despite the 
irregular grid. 



#v 


^^max 

err^2 

err^i 

eiTL°- 

^min 

#Newton 

nonlinear scheme 

324 

2.E-04 

1 

3.99E-02 

1 0.404 

1.42E-02 

8.92E-04 

1148 

linear scheme 

324 

2.E-04 

1 

1 3.47E-02 

0.377 

2.01E-02 

-1.49E-02 

259 


Table 3. Kershaw grid. Nonlinear and linear scheme, with an 
anisotropic tensor. 


Denoting by rc = tt exp {g{x — 5 )) the long-time asymptotic of u defined by (110), 
then the relative entropy of a function u : > K-|- w.r.t. w is defined by 

(111) E^{u)= J l^ulog ^ — u-I-da;. 

It is simple to verify that 

Therefore the decay of the free energy is equivalent to the decay of the relative 
entropy. Note that A™ is undehned (or is set to -boo) if u < 0 on a positive 
measure set. It is well known (see e.g. [36, 77]) that the relative entropy 
converges exponentially fast towards 0 as t tends to -boo. Exponential convergence 
results in the discrete setting were proved for instance in [39, 38, 17] in the case of 
a monotone discretization of dissipative equation (see also [16]). In order to check 
this asymptotic behavior at the discrete level, we introduce the discrete relative 
entropy E^{u) defined for all nonnegative u = (uk,'Us)„s € W-d (i.e., such that 
UjS >0 for all /3 G A1 U V) by 

(112) E^{u)= ^ mpi^plog 

/3GAIUV 

The exponential convergence towards equilibrium is recovered as it appears clearly 
on Figure 5. 


Up 


’{xp) 


-Up+ w(xp)"j . 
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Figure 5. Left: Kershaw mesh. Right: Evolution of the relative 
entropy 1 1 —>■ t)) on a logarithmic scale in function of time. 


5.2.2. Test 2: Porous medium equation with Dirichlet boundary condition. In this 
section, we apply our scheme to the case of the anisotropic porous medium equation 

(113) dfU — V- (A V(it^) ) = 0 in 

for different choices of functions ry and p with r]{a)p'{a)da = namely 

(a) 77 ( 11 ) = 2u^ and p(u) = log(u), 

(b) 77 ( 77 ) = 2u and p{u) = u, 

(c) 77 ( 77 ) = 1 and 77 ( 77 ) = 77 ^. 

For the choice (a), the function rj is strictly convex. Therefore, the rigorous gra¬ 
dient flow structure of the problem corresponding to this choice of mobility function 
77 is unclear [43]. The pressure function p is singular near 0, hence Lemma 3.7 im¬ 
plies that the corresponding scheme is positivity preserving. 

The choice (b) with a linear mobility corresponds to the now classical setting 
highlighted in [86, 77]. 

Finally, the choice (c) corresponds to the usual approach for discretizing the 
porous medium equation. The corresponding scheme enters into the framework 
of [48], where its convergence is proved. 

The problem is closed here with Dirichlet boundary conditions (destroying by 
the way the gradient flow structure but not the convergence of the scheme). 

Comparison with a one-dimensional analytical solution. The numerical convergence 
of the three schemes has first been compared thanks to the following analytical 
solution (again built in 1-dimension), 

(114) u{x,y,t) = max{2Ct - X , 0), y{{x,y),t) & Qti- 

Note that (114) is the unique weak solution corresponding to the initial condition 
uo(x, y) = u{x, 77 , 0) and to the Dirichlet boundary condition ud{x, y, t) = u{x, y, t) 
on dTl X (0,tf). Our numerical convergence study makes use of the family of trian¬ 
gular meshes already used for Test 1. Once again, the final time t{ is fixed to 0.25 
and an anisotropic tensor is given by /a, = 1 and ly = 10. 

We observe in Tables 4-6 that second order convergence is destroyed for all the 
three schemes because of the lack of regularity of the exact solution. As expected, 
the discrete solution corresponding to the choice (a) remains positive while the 
discrete solutions to the schemes corresponding to the choices (b) and (c) suffer of 
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h 

#v 

^^init 

^^max 

err £,2 

rate 

err_£i 

rate 

err/^oc 

rate 

-fimin 

T^Newton 

0.306 

37 

0.001 

0.01024 

0.523E-02 


0.997E-03 

- 

0.105E+00 


0.000 

479 

0.153 

129 

0.00025 

0.00256 

0.205E-02 

1.352 

0.344E-03 

1.535 

0.522E-01 

1.013 

0.000 

1143 

0.077 

481 

0.00006 

0.00064 

0.898E-03 

1.190 

0.123E-03 

1.490 

0.259E-01 

1.012 

0.000 

2218 

0.038 

1857 

0.00002 

0.00016 

0.380E-03 

1.240 

0.417E-04 

1.554 

0.128E-01 

1.012 

0.000 

5652 


Table 4. Test 2: Choice (a) of mobility and pressure functions, 
convergence towards (114). 


h 

#v 

Atinit 

Aimax 

err £,2 

rate 

err£,i 

rate 

err£,c« 

rate 

Umin 

^Newton 

0.306 

37 

0.001 

0.01024 

0.769E-02 

- 

0 . 210 E -02 


0.645E-01 


-0.032 

138 

0.153 

129 

0.00025 

0.00256 

0.263E-02 

1.546 

0.613E-03 

1.775 

0.326E-01 

0.983 

-0.017 

383 

0.077 

481 

0.00006 

0.00064 

0.897E-03 

1.554 

0.173E-03 

1.823 

0.164E-01 

0.996 

-0.009 

1246 

0.038 

1857 

0.00002 

0.00016 

0.306E-03 

1.551 

0.481E-04 

1.849 

0.821E-02 

0.996 

-0.005 

4234 


Table 5. Test 2: Choice (b) of mobility and pressure functions, 


convergence towards (114). 


h 

#v 


^^max 

err £^2 

rate 

err£,i 

rate 

err£,«> 

rate 

Umin 

^Newton 

0.306 

37 

0.001 

0.01024 

0 . 116 E -01 

- 

0.371E-02 


0.764E-01 


-0.065 

148 

0.153 

129 

0.00025 

0.00256 

0.423E-02 

1.461 

0.116E-02 

1.672 

0.388E-01 

0.977 

-0.039 

436 

0.077 

481 

0.00006 

0.00064 

0.149E-02 

1.501 

0.337E-03 

1.788 

0.233E-01 

0.737 

-0.021 

1438 

0.038 

1857 

0.00002 

0.00016 

0.524E-03 

1.513 

0.932E-04 

1.856 

0.129E-01 

0.856 

-0.010 

4912 


Table 6. Test 2: Choice (c) of mobility and pressure functions, 
convergence towards (114). 


undershoots. The choice (b) appears to be both cheaper and more accurate than 
the choice (c), and the amplitude of the undershoots is smaller. 



Figure 6. Test 2. Coarsest grid. Discrete unknown (us)sgv„ and 
its iso-values. Choice (a) (left) and (c) (right) for 77 and p. 


Figure 6 illustrates the iso-values of the piecewise affine functions defined on 
the triangular mesh Ai reconstructed thanks to its nodal values (u")ggy for the 
coarsest triangle grid at the final time tf. For the choice (a) of the mobility and the 
pressure (left), the iso-values are chosen from 0 to 0.025 by step of 0.05 and then 
from 0.1 to 0.5 by step of 0.1. For the choice (c) of the mobility and the pressure 
(right), the iso-values are taken from —0.025 to 0.025 by step of 0.05 and also from 
0.1 to 0.5 by step of 0.1. 





































38 


CLEMENT GANCES AND CINDY GUICHARD 


Comparison with a two-dimensional analytical solution. We test our approach on 
the two-dimensional analytical solution 


(115) 


a(a; - 0.5)2-h/3(j/- 0.5)2 


in Qtt 


of the anisotropic porous medium equation (113), where tf has been set to 0.25. 
The permeability tensor is still assumed to be diagonal with l^ = 0.1 and ly = 10, 
and we have set a = and /3 = The problem is closed with Dirichlet 

boundary conditions and the initial condition corresponding to (115). The results 
are gathered in Tables 7-9. 


h 

#v 

Atinit 

Atmax 

err^Q 

rate 

err 2^1 

rate 

err^^ 

rate 

Umin 

#Newton 

0.306 

37 

0.001 

0.01024 

0.270E-02 


0.114E-02 


0.120E-01 


0.000 

89 

0.153 

129 

0.00025 

0.00256 

0.942E-03 

1.517 

0.381E-03 

1.578 

0.473E-02 

1.341 

0.000 

223 

0.077 

481 

0.00006 

0.00064 

0.293E-03 

1.688 

0.119E-03 

1.686 

0.163E-02 

1.534 

0.000 

805 

0.038 

1857 

0.00002 

0.00016 

0.802E-04 

1.868 

0.334E-04 

1.828 

0.461E-03 

1.826 

0.000 

3142 


Table 7. Test 2: Choice (a) of mobility and pressure functions. 


convergence towards (115). 


h 

#v 

Atjnit 

Atniax 

err £^2 

rate 

err^^i 

rate 

err/, CO 

rate 

'i^miii 

T^Newton 

0.306 

37 

0.001 

0.01024 

0.645E-02 

- 

0.271E-02 


0.271E-01 


-0.027 

99 

0.153 

129 

0.00025 

0.00256 

0.202E-02 

1.676 

0.828E-03 

1.709 

0.992E-02 

1.447 

-0.008 

237 

0.077 

481 

0.00006 

0.00064 

0.604E-03 

1.742 

0.246E-03 

1.753 

0.341E-02 

1.540 

-0.002 

801 

0.038 

1857 

0.00002 

0.00016 

0.161E-03 

1.905 

0.667E-04 

1.882 

0.966E-03 

1.821 

0.000 

3140 


Table 8. Test 2: Choice (b) of mobility and pressure functions. 


convergence towards (115). 


h 

#v 

Atjnit 

Atjuax 

err/,2 

rate 

err /,1 

rate 

err/, CO 

rate 

'i^min 

^i^Newtoii 

0.306 

37 

0.001 

0.01024 

0.102E-01 

- 

0.448E-02 


0.575E-01 


-0.046 

128 

0.153 

129 

0.00025 

0.00256 

0.321E-02 

1.661 

0.134E-02 

1.739 

0.194E-01 

1.569 

-0.011 

250 

0.077 

481 

0.00006 

0.00064 

0.933E-03 

1.783 

0.383E-03 

1.811 

0.553E-02 

1.808 

-0.003 

810 

0.038 

1857 

0.00002 

0.00016 

0.244E-03 

1.933 

O.lOlE-03 

1.927 

0.147E-02 

1.914 

-0.001 

3140 


Table 9. Test 2. Choice (c) of mobility and pressure functions. 


convergence towards (115). 


As expected, the choice (a) leads to a positivity preserving scheme, contrarily to 
the choices (b) and (c). Moreover, the scheme (a) is the most accurate and does 
not come with an additional cost. 


5.2.3. Test 3: Porous medium equation with drift. In this third test case, we have 
set rj{u) = M on IR+ and p(u) = 2u and 5 = 1, leading to the degenerate problem 

(116) i9tu - V-(A (V(u2) - Mg) ) = 0 inQtj. 

The problem is endowed with Dirichlet boundary conditions. The tensor A is chosen 
to be diagonal with lx = ^ and ly = 100. We compare the results obtained by (43) 
with those obtained using, instead of (43c), this particular definition of the fluxes 


(117) F,,s(m”) = 


E' 

s'GV„ 


«,,(«)2-(Mr0") + 


■E' 

s'G Vk 




Vk e M, Vs e Vk- 
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The resulting scheme is called the quasilinear scheme. The numerical convergence 
of both schemes has been compared on the sequence of triangular meshes already 
used in the previous tests, thanks to the following analytical solution (again built 
in 1-dimension), 

(118) u{x,y,t) = m.ax{f3t — X , 0), V((a;,y),t) S x (0,tf), 

with /3 = -I- g). The profile (118) is the unique weak solution corresponding to 

the initial condition uq{x, y) = u(x, y, 0) in fl and the Dirichlet boundary condition 
UD{x,y,t) = u(x,y,t) on dO, x (0,tf). 


h 

#v 

^^init 

^^max 

err £,2 

rate 

err£,i 

rate 

errL~ 

rate 

Umin 

^Newton 

0.306 

37 

0.001 

0.01024 

0.130E-01 

- 

0.423E-02 


0.890E-01 


-0.046 

187 

0.153 

129 

0.00025 

0.00256 

0.495E-02 

1.398 

0.133E-02 

1.675 

0.496E-01 

0.843 

-0.032 

552 

0.077 

481 

0.00006 

0.00064 

0.184E-02 

1.428 

0.397E-03 

1.741 

0.283E-01 

0.808 

-0.017 

1609 

0.038 

1857 

0.00002 

0.00016 

0.660E-03 

1.479 

0.116E-03 

1.771 

0.145E-01 

0.970 

-0.009 

5586 


Table 10. Test 3. Nonlinear scheme (43). 


h 

#v 

^^init 

^^max 

err £,2 

rate 

err£,i 

rate 

err£,co 

rate 

Umin 

^Newton 

0.306 

37 

0.001 

0.01024 

0.154E-01 

- 

0.568E-02 


0.939E-01 


-0.068 

193 

0.153 

129 

0.00025 

0.00256 

0.671E-02 

1.201 

0.213E-02 

1.416 

0.613E-01 

0.615 

-0.048 

642 

0.077 

481 

0.00006 

0.00064 

0.271E-02 

1.309 

0.702E-03 

1.600 

0.326E-01 

0.910 

-0.027 

2178 

0.038 

1857 

0.00002 

0.00016 

0.104E-02 

1.384 

0.212E-03 

1.725 

0.170E-01 

0.938 

-0.015 

7365 


Table 11. Test 3. Quasilinear scheme, fluxes defined by (117). 


Here again, the convergence orders of both scheme are similar, but strictly lower 
than 2 because of the lack of regularity of the exact solution. Both schemes violate 
the positivity of the solution in this case, but the amplitude of the undershoots is 
smaller for the nonlinear scheme. There is no contradiction here with Lemma 3.7 
since p is not singular at u = 0. Our nonlinear scheme is slightly more accurate, 
produces undershoots with a smaller amplitude, and is cheaper than the quasilinear 
one. 


5.2.4. Test 4- ^ heterogeneous test case. The last test aims to illustrate the ability 
of the scheme to deal with heterogeneous situations. Motivated by an application to 
complex flows in porous media (see for instance [35, 32]), we test the nonlinear VAG 
scheme in a slightly more complicated configuration where both the permeability 
tensor A and the pressure function p depend on a: in a discontinuous way. More 
precisely, the domain fl = (0,1)^ is made of two open subdomains Hi (the drain) 
and H 2 (the barrier) with H = Hi U H 2 and Hi n H 2 =0 (see Figure 7 for a 
representation of Hi and H 2 ). The permeability tensor and the pressure function 
are defined by 


and 


A(a:) 


p{u,x) 


Ai =Id 


if a; S Hi, 

(1 



A 2 = 

if a; e H 2 , 

VO 

0.01 J 


ipiiu) = 

31og(u) 

if a; G Hi 

\P2(W) = 

log(M) 

if a; G H 2 
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The mobility function is linear and does not depend on x, i.e., r]{u) = u. For the 
sake of simplicity, we have set 1^ = 0. At the interface J between fli and ^21 
the flux and the pressure are assumed to be continuous, i.e., denoting by Ui the 
restriction of u to fli and by rii the normal to J outward w.r.t. fli, we require 

(119) 3AiV mi • ni + A 2 VM 2 • 412 = 0, and Pi{ui) = P 2 {u 2 ) onJx{0,tf). 


The problem is complemented with the boundary conditions 

• p{u, a:) = 0 (hence u = 1) on the bottom boundary, 

• p{u, x) = —4 (hence u ~ 0.018) on the top boundary, 

• AVtt • n = 0 on the lateral boundaries. 

The initial data is chosen at equilibrium, with p(uo,x) = —4 in the whole Q. 
Existence and uniqueness for this problem follow from the analysis carried out 
in [28]. 

Since u is discontinuous across J (in opposition to the pressure p{u, ■) follow¬ 
ing (119)), it is natural to choose p rather than u as the primary variable of the 
numerical scheme (cf. [62, 55], we refer to [23] for an alternate strategy that improves 
robustness) in order to avoid the complex treatment of the jump condition (119) 
at the interface performed for instance in [28, 29, 47, 24]. 

The mesh A4 is assumed to be compatible with the geometry of fl, in the sense 
that K S A4 is either contained in or 112, but J n k = 0 for all k S A4 (cf. 
Figure 7). Define the functions : M —>■ (0, 00 ) as the inverse of p{-,x^) for all 
K G Ai. The subset of V made of vertices belonging to the top or bottom boundaries 
where Dirichlet boundary conditions hold is denoted by Vext- We also make use of 
the notations Vint = k’ \ Vext and V^.int = Vint n V^ for k G A4. The scheme (43) 
expressed with p as a primary variable consists in finding p = (p]l,p")^ s n W'p, 2 ^t 
such that for all n > 1, 




sev,. 


V- u^{Ps)-u^{Ps 
> . -TT- 


At 

kGA4s 

^k”s+^s”.=0, 




kGA4b 




= 


s'GV,, 
W«(p;:) +M«(p”) 


Vk G A4, 

0, Vs G Vint, 

Vk G A4,Vs G VK,int, 
Vk G A4, Vs G Vk, 

Vk G A4, Vs G Vk. 


We observe on Figure 8 that the results on the triangular mesh (with 481 nodes) 
and on the cartesian mesh (with 289 nodes) are similar. Moreover, the numbers of 
Newton-Raphson iteration needed to compute both solutions are of the same order. 


Appendix A. Some lemmas related to the VAG discretization 

This appendix gathers lemmas on some properties of the VAG discretization 
that are independent of the continuous problem (and thus of the scheme). In 
what follows, V = {A4,T) denote a discretization of D as prescribed in §2.1.1, and 
TTp, 7 r;\ 4 , TTx) and V 7 - are the corresponding reconstruction operators. 
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Figure 7. Test 4. Illustration of the two sub-domains: the drain 
rii in blue and the harriers ^2 in red. Left: cartesian grid. Right: 
unstructured triangular grid. 


Lemma A.l. For k € M., let be the matrix defined by (34), 

then there exists C depending only on A, Oj- and £x> (but not on k) such that 
Cond2(A«) < C. 

Proof. Following [26, Lemma 3.2], there exist Ci,C 2 > 0 depending only on Oj- and 
tj) such that, for all u G W-d and all k G Ai, one has 

in63/S(/^) V—^ , s 2 i\^m, ii2 ^ m6cLS(K) v—> , ,2 

(h z2 < W^tuWl^^) < C'a .2 z2 ’ 

^ sGV,, ' sGV,^ 

where h^ denotes the diameter of the cell k G Ai. As a consequence, one has 


A*Ci 


meas(«:) 


|dKLt|^ < • Ak^kM = / A'V'ru-'Vrudx < X*C2 ^(---v(^\Sk.u\ 


{K? 

Since the application : IFp 
meas(tt;). .2 


meas(«:) 

7 ^' 


is onto, we deduce that 


A*Ci- 


IV . , * ^ meas(K), ,, 

vY <v A^v < X*C 2 AM t'P 


Vu € 


{h^r ' ' - " - " (h.y 

and thus that Cond2(AK) < ^ . □ 

Lemma A.2. There exists C > 1 depending only on A, Oj- and tx> such that, for 
all K G At and all v = (us)sgv,, ^ 


<Cv- A^v. 

sGV,, Vs'GV,, / 

Proof. Denoting by || • ||q the usual matrix 5 -norm, one has 

E ( E K^'i) < iiA^iiiiiir 

sGV,. Vs'GV,, / 

Since the dimension of the space is bounded by £t>, there exists Ci depending 
only on €-p such that ||Ak||i < Ci||A,.;|| 2 , so that 
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Figure 8. Approximation of u{p,x) at times t = 0.05, t = 0.2, 
and t = 1 for the two different meshes. 


On the other hand, since is symmetric definite and positive, one has 


V ■ A^v > 


Cond2(A«) 


Using Lemma A.l, we obtain that there exists C 2 > 0 depending only on A, Oj- 
and Ixi such that 


(121) 


V ■ A^v > C'2|jA«,||2|t>p. 
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Putting (120) and (121) together, we conclude the proof of Lemma A.2 by choosing 

c=a. □ 


Lemma A.3. Let k G M and A,^ = {a'^g,)s^s'ev^ £ 




be the matrix defined 


by (34). Let = (/i^.^) 


sGV,^ 


and V G Wxi, then 


sGV,, s'GV„ ^ " sGV„ Vs'GV,, / 


Proof. Using ab < ^ we obtain that 


X! (^s “ i’K)t^K.,sas,s'h-'ipivs' - v^) 

< max(^«,s)^ X! - T^||a,>||Ws/ - 


sGV,, s'GV,, 


< 


sGV,, 

sGV,^ s'GV,^ 

maxsgv,,(MK,s)^ ( ST^ I K I 1 / n2 

sGV„ Vs'GV,, / 

maxsgv,,(AiK,s)^ i « i 1 c ' ^2 

+- 2 -^^ ^ • 

s'GV,, \sGV„ / 

One concludes the proof of Lemma A.3 by noticing that, since A„ is symmetric, 
the two terms in the right-hand side of the above inequality are equal. □ 


Lemma A.4. There exists C depending only on 9j- and (.x) sueh that 

meas(T) < meas(K) < C'meas(T), Vk G Ai, VT G Twith T C n. 

Proof. Let k G Ai, then there exist Ti,..., simplexes, with r = \i d = 2 and 
r = 2#£k if d = 3, such that 

T'k. 

[Jt, = k, t, n t, = 0 if i ^ j. 

i=l 

The Euler-Descartes theorem ensures that r < 4(£'p — 1) if d = 3. 

If Ti and Tj share a common edge, one gets that 

meas(ri) < 9‘^ meas(Tj). 

Let io,ii £ {!)• ■•)*’«} be arbitrary but different, we deduce from the previous 
inequality the following non-optimal estimate: 

meas(Ti(,) < meas(TiJ. 

Let fniax be such that meas(Ti^^^) = maxi<i<r meas(Ti), then 

meas(N) < r measiTi^^^) < A{(.x — l)d^^^'°~^^'^meas(Ti), Mi G {1,..., r}. 

□ 


We state now a slight generalization of [26, Lemma 3.4], where the same result 
is proven in the particular case q = 2. The straightforward adaptation of the proof 
given in [26] to the case g ^ 2 is left to the reader. 
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Lemma A.5. There exists C depending only on I-d and Oj- defined in (24) and (21) 
respectively such that, for all v G VLp and all q G one has 

hw - + hvv - TTMv\\L,(n) ^ C'VllVr^’l!L,(n)- 


Lemma A.6. Let D be a discretization of Q. as introduced in such that 

Qd > 0, then there exist Ci > 0 depending only on q, 9j- and £x> and C 2 depending 
moreover on Qd such that 


(122) C'i||7rx)v||L<!(n) < lkri’l!L‘!(n) < C'2||7r-DL’||L<!(n), yvGW-D- 


Proof. Let T be a reference tetrahedron, and let : T —)■ K be an afhne function 
with nodal values Vi, i G {1,.. .A}, then for all q > 0, there exists C depending on 
q such that 


1 

C 


E 




9 

L9(f) 


i=l 


Therefore, using classical properties of the affine change of variable between sim- 
plexes, one gets the existence of C depending only on q, Op, and f-p such that, for 
all V G Wt>, 


(123) ^ ineas(K) ( |?;«|« + ^ 


k<^M 


sGV,^ 


^ ^ C' E meas(K) ^ j • 

V sGVk / 


On the other hand, one has 

kGA4 sGV 

A classical geometrical property and (29) yield 

f d 

(124) Wk < meas(K) =d TTpe,^{x)dx < —Vk G M, 

Jn Qv 

and similarly 

/» ^ 

mg < d 7r7-es(£c)da: < —m^, Vs G V. 

Jn Cx> 

Notice now that the following geometrical identity holds: 

d / 7r7-es(a;)dx = ^ meas(T), Vs G V. 

t/ O T ^ 7 ~ 

x^GdT 

Lemma A.4 yields the existence of C > 0 depending on 9j- and £■£> such that 
^ ^ meas(«:) < d j TT'fes{x)dx < ^ meas(/c), Vs G V, 


kGMs 


kGMs 


and the result of Lemma A.6 follows. 


□ 
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Lemma A.7. Let D be a discretization of Ll as introduced in such that 

Qd > 0, then, for all q G one has 




Ct> J 


1/9 


\TTx,V\ 


Li(Q) 


Vv G W-D- 


Proof. Let v = ggy G tLp, then it follows from (124) that 

Ik^^llL(n) = meas(K) 


kGM 


^ {£) E ^ (^) ■ 


□ 


Lemma A.8. Let v = {vk.,Vs)^ ^ G Wt> be such that > 0 for all fi G M UV, 
and define v = (U„,Ws)k s ^ 


Vs = 0 , 


Vk. = max Vk., max Vs' , 
s'GV^ 


Vs G V,Vk G M. 


Then there exists C depending only on Op, (.x) o.nd Qx such that 

hMv\\Li{n) < C'IKx>ii|lLi(n) ■ 

Proof. Let v G Wx be a vector with positives coordinates, and let v be constructed 
as above. It follows from the construction of v that 


Vk < Vk, + '^ Vg, Vk G M, 


sGVk 


whence, applying (123) with q = 1, one gets 

IKaiHIIlhsi) — C'|| 7r7-L’||Li(n). 
The result now directly follows from Lemma A.6. 


□ 


Lemma A.9. Let u = (u^,M s)kgA 4 .sgv G Wx, then for all k G M., we define 
du = G Wx by 


SgU = 0 and SkU = max \uk — ' 

s'GVk 


Vk G M, Vs G V, 


then, for all q G [l,oo], there exists C depending only on q. Op, and ix such that 
(125) < Chr\\VTu\\Li{n)- 

Proof. Let k G A4 and s G Vk, then there exists a simplicial sub-element T G T oi 
K G Ai such that Xk and Xg are vertices of T. Then it follows from classical finite 
element arguments (see e.g. [40, 46]) that 

meas(r)^/«|uK - Ms| < c^EL|j < Chrll VrM|lL‘!(K), 

Pt 

where c depends only on the dimension d and on q, while C depends additionally 
on 6p. Thanks to Lemma A. 4, we get the existence of C depending on d, q, Op and 
^x such that, 

meas(K)^'^'^|uK - Usl < Chp |1 VrM|li,(K), Vk G M, Vs G V^. 

Summing over k G M. provides that (125) holds. □ 


Vk G M, Vs G V^. 
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